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ABSTRACT 

We implement physically motivated recipes for partitioning cold gas into different 
phases (atomic, molecular, and ionized) in galaxies within semi-analytic models 
of galaxy formation based on cosmological merger trees. We then model the 
conversion of molecular gas into stars using empirical recipes motivated by recent 
observations. We explore the impact of these new recipes on the evolution of 
fundamental galaxy properties such as stellar mass, star formation rate (SFR), 
and gas and stellar phase metallicity. We present predictions for stellar mass 
functions, stellar mass vs. SFR relations, and cold gas phase and stellar mass- 
metallicity relations for our fiducial models, from redshift 2 ~ 6 to the present 
day. In addition we present predictions for the global SFR, mass assembly history, 
and cosmic enrichment history. We find that the predicted stellar properties 
of galaxies (stellar mass, SFR, metallicity) are remarkably insensitive to the 
details of the recipes used for partitioning gas into Hi and H 2 . We see significant 
sensitivity to the recipes for H 2 formation only in very low mass halos, which 
host galaxies that are not detectable with current observational facilities except 
very nearby. The properties of low-mass galaxies are also quite insensitive to 
the details of the recipe used for converting H 2 into stars, while the formation 
epoch of massive galaxies does depend on this significantly. We argue that this 
behavior can be interpreted within the framework of a simple equilibrium model 
for galaxy evolution, in which the conversion of cold gas into stars is balanced on 
average by inflows and outflows. Star formation in low mass galaxies is strongly 
self-regulated by powerful stellar driven outflows, so the overall galaxy-scale star 
formation efficiency is nearly independent of the H 2 depletion time. Massive 
galaxies at high redshift have not yet had time to come into equilibrium, so the 
star formation efficiency is strongly affected by the H 2 depletion time. 
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1 INTRODUCTION 


While the ACDM (Cold Dark Matter plu s cosmological 
constant A) model ( Blumenthal et aUl984l ) now provides 
us with a well-motivated framework for predicting the 
abundances and properties of dark matter halos and the 
large scale structures in which they are embedded, all 
galactic or larger scale simulations must rely on “sub¬ 
grid” recipes in order to treat processes such as star 
formation and stellar feedback. Cosmological simulations 
are unable to directly resolve individual stars or, usually, 
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even Giant Molecular Clouds (CMC). In order to model 
the conversion of cold gas into stars, up until recently, 
both numerical and semi-analytic cosmological simula¬ 
tions typically utilized a very simple empirical sub-grid 
recipe based on observations most famou sly by ISchmidtl 
(ll959l . Il963h and iKeririiciiul (Il989l . 1 19981 ) (often referred 
to as the “Kennicutt-Schmidt” (KS) relation). These ob¬ 
servations showed that the surface density of star forma¬ 
tion Esfr was proportional to the surface density of cold 
gas to a power Nks- Observations also showed that the 
efficiency of star forma tion dropped rapidl y b elow a criti¬ 
cal gas surface density ((Martin fe Kennicuttl200ll ). There 
has been debate about whether this critical surface den- 
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sity is best desc ribed in terms of a Toomre stability cri¬ 
terion dToomrdll964h or a constant critical density, and 
indeed about the physical origi n of this critical density 
JSchavell2004l : ILerov et, al.ll2008l 'l. 

From the pioneering work of iKatd (1 19921 1 up until 
recently, cosmological simulations of galaxy formation, 
both numerical and semi-analytic, have implemented a 
star formation recipe in which “cold” gas (typically with 
T ^ 10 4 K) with volume density p gas is assumed to form 
stars at a rate per unit volume: 

P* ^*Pgas (1) 

with N ~ 1.5 and e* usually treated as a free param¬ 
eter, tuned to match the observed Kennicutt relation. 
A common variant assumes p* oc p gas /tfi, which is ap¬ 
proximately equivalent because the local free-fall time 
tfr oc p~ 0 ' 5 . Motivated by the observational evidence de¬ 
scribed above, many modelers incorporated either a crit¬ 
ical surface density or volume density into their star for¬ 
mation recipe, which proved to be important in order to 
reproduce the observed high gas fractions in low-mass 
galaxies. 

Beginning about a decade ago, our understanding 
of how star formation on ~ 100 pc-kpc scales de¬ 
pends on lo cal con ditions began to undergo a revolu¬ 
tion. IWong fe Blit j d2002i ') showed that the correlation 
between Esfr and the surface density of molecular gas 
Tjh 2 was stronger than that between Esfr and the total 
gas density E gaa in molecule rich galaxies. In the past 
five years, this field has advanced rapidly with the avail¬ 
ability of galaxy-wide, high resolution maps of the star 
formation and multi-phase (Hi and H 2 ) gas in reasonably 
large samples of nearby galaxie s, e.g. from the TH INGS 
(The HI nearby galaxy survey: Walter et al.ll2008r ) com¬ 
bined with BIMA SONG (BIMA survey of Nearby Galax¬ 
ies; I Heifer et alJ 120039 an d HERACLES (H ERA CO- 
Line Extragalactic Survey; ILerov et alJ 120031 . Bas ed o n 
these observations, it has been show n jBigieLgLaL [20081; 
ILerov et al .1120081 : iBigiel et, al.ll201 ll : ISchruba et al. 1201 If ) 
that, when averaged over scales of ~ 700 pc, the star 
formation density is tightly correlated with the molec¬ 
ular gas density to a nearly linear power, and that 
there is almost no correlation between Esfr and the 
density of atomic gas, so that the correlation between 
Esfr and E gaa (the traditional KS relation) breaks down 
badly in the Hi-dominated parts of galaxies (typically in 
galaxy outskirts). These results highlight the importance 
of modeling the partition of gas into different phases, 
i.e. atomic vs. molecular, which has not been attempted 
in most cosmological simulations of galaxy formation to 
date. 

At the same time, there has been significant progress 
in understanding and modeling the formation of molec- 
ular hydrogen and star f ormat ion on galactic scales. 
iBlitz fe Rosolowskvl d2004l . 120069 showed that the frac¬ 
tion of atomic to molecular gas in a sample of nearby disk 
galaxies was tightly correlated with the midplane pres¬ 
sure (determined by the density of both stars and gas), 
and this resul t has b een confirm e d in larger samples such 
as THINGS (ILerov et al.l 120089 . iRobertson fe Kravtsov! 
(|20089 implemented low-temperature (T < 10 4 K) cool¬ 
ing, photo-dissociation of H 2 , and an H 2 -based SF recipe 


in hydrodynamic si mulations of iso l ated di sk galaxies 
of various masses. iKrumholz et al.1 J2009b9 presented 
analytic models for the formation of H 2 as a func¬ 
tion of total gas density and metallicity, supported 
by numerical simulat i ons w ith simplified geometries 
(IKrumholz et al.1120081 . l2009 ah emphasizing the impor¬ 
tance of metallicity as a controlling para meter in H 2 
formation. iGnedin fe Kravtsov! d20ld . 1 20 Ilf ) included de¬ 
tailed chemistry and low temperature cooling as well 
as a simplified treatment of radiative transfer and an 
H 2 -based SF recipe in cosmological “zoom-in” Adaptive 
Mesh Refinement (AMR) simulations of small regions, 
and presented analytic fitting functions to their results 
as a function of total gas density, m etallicity, and the 
streng th of the local UV background. IChristensen et ahl 
(120129 used a similar approach to implement chemistry 
and simplified radiative transfer in Smoothed Particle 
Hydrodynamics (SPH) zoom-in simulations of galaxy¬ 
sized regions, which include a blast-wave treatment of 
supernova feedback. 

A somewhat dif ferent view has been presented by 
lOstriker et al.l (120109 . who propose that heating of the 
Interstellar Medium (ISM) by the stellar UV background 
plays a key role in regulating star formation. In their 
model, the thermal pressure in the diffuse ISM, which 
is proportional to the UV heating rate, adjusts until it 
balances the midplane pressure set by the vertical grav¬ 
itational potential. This could provide an explanation 
for the strong empirical correlation b etween H 2 fraction 
and d isk midplane pressure found bv lBIitz fe Rosolowskvl 

(I2006IL 

Although detailed simulations are crucial in order 
to understand the complex physical processes involved, 
extremely high resoluti on is required in order t o obtain 
reliable results (see e.g. IChristensen et al.ll2012r) . imply¬ 
ing that it will be feasible to simulate only small num¬ 
bers of galaxies with these techniques for the next few 
years. Meanwhile, large surveys of cold gas in nearby and 
distant galaxies with new and upcoming facilities such 
as the Atacama Large Millimeter/submillimeter Array 
(ALMA) and the SKA (Square Kilometer Array) and its 
pathfinders are already being planned and pilot projects 
are underway. As a result, it is important to develop 
computationally efficient techniques that can incorporate 
physically motivated treatments of gas partioning into its 
atomic, molecular, and possibly ionized phases and H 2 - 
based star formation recipes into simulations of cosmo¬ 
logical volumes. 

Semi-analytic models (SAMs) provide an alterna¬ 
tive approach to this problem. In semi-analytic merger 
tree models, a merger tree represents the formation and 
growth of a dark matter halo that is identified at some 
redshift of interest; these merger trees may be extracted 
from dissipationless N-body simulations o r created us¬ 
ing analytic techniqu es (e.g. ISomerville fe Kolattl 1 19991 : 
Parkinson et al. [20089 . Simplified but physically moti¬ 
vated recipes are used to track the rate of gas cool¬ 
ing into galaxies, and these recipes have been tested 
against fully numerical hydrodynamic simulations (e.g. 
iHirschmann et al.ll20 12aL These models use angular mo¬ 
mentum based arguments to track the radial sizes o f 
forming disks (IMo et al.l 1 19981 : ISomerville et al.l l2008b9 , 
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and can then implement recipes for how cold gas is con¬ 
verted into stars, and how energy and momentum from 
massive stars and supernovae is returned to the Interstel¬ 
lar Medium (ISM). This “feedback” from stars and SNae 
is assumed to drive large-scale winds that can remove gas 
from the galaxy. The production, ejection, and recycling 
of metals is also tracked. Thus our existing semi-analytic 
modeling framework provides the main quantities (total 
gas density in disks, gas metallicity) needed to implement 
physically motivated recipes for partitioning gas into an 
atomic and molecular component and then implementing 
an H 2 -based star formation recipe. 

Several efforts along these lines_ h ave already been 
presen ted in the literature. lObreschkow fe Rawlingsl 
( 2009 1 implemented a prescription to estimate the H 2 
fr action based on the empirical pressure-based recipe 
of lBIitz fe Rosolowskvl (hood ) applied in post-processin g 
to th e Millennium simulations of IDe Lucia fc Blaizotl 
(|2007l 'l. However, in this approach, the star formation 
in the simulations was still based on a traditional KS 
recipe using the total gas density, not sel f -consi s tently 
on the estimated H 2 gas density. IfyT et all (<2010l . 120121 1 
modeled the partitioning of gas into Hi and H 2 in ra¬ 
dial bins in each ga laxy, using both the m etallicity- 
dependent recipes of K rumholz et a l. (12009b: ) an d the 
pressure-based recipe of lBIitz fe Rosolowskvl d2006h . and 
self-consistently implemented an H 2 -based star forma¬ 
tion reci pe, within the semi-analytic modeling frame¬ 
work of Iguo et all ll201lll . iLaeos et al.l (l2011bl lah also 
estimated gas partitioning into an atomic and molec¬ 
ular component, and implemented an H 2 -based star 
formatio n recipe, within the GALFORM sem i-analytic 
models (IPaugh et alJ I 2 OO 5 I : iBower et al.l l200fil l. Similar 
modeling efforts utilizing a somewhat more simplified 
framework (i.e., only the mass accretion history of the 
largest progenitor is tracked, rat her th an the full merger 


the pressure-based Blitz & Rosolowskv ( 

20060 approach, 

and Krumholz & Dckel 

2012h using the 

Krumholz et al. 


It is already clear that the results of this kind 
of exercise may depend on the other ingredients of 
the modeling, in particular on the treatment of stellar 
feedback, chemical evolution, and potentially on feed¬ 
back from Active Galactic Nuclei (AGN). In this work, 
we present new models that incorporate a metallic¬ 
ity or pressure oriented treatment of atomic-molecular 
gas partitioning and an H 2 -based star formation recipe 
within the semi-analytic modeling framewor k developed 
by the Santa Cruz group (ISomerville fe Primackl 1 19991 : 

I Somerville et al.ll200ll. l2008al . 1201 

The current generation of semi-analytic models (in¬ 
corporating some form of “quenching” in massive ha¬ 
los, e.g. from AGN feedback) has been fairly success¬ 
ful at reproducing a variety of galaxy observations, 
but suffer from generic problems. Both the successes 
and problems seem to be common to the semi-analytic 
models developed by many different groups as well as 
to lar ge-volume cosmological hydrodynamic simulations 
(see ISomerville fe Da ve 2014, for a discussion). Signifi¬ 
cant successes include the ability to match the observed 
stellar mass function or luminosity functions from the 


UV to the NIR at z = 0, while simultaneously match¬ 
ing the gas fraction as a function of stellar mass for 
nearby disk ga laxies (e.g. ISomerville et al.ll2008al . 120121 : 
iLu et al.l 120141 ). Observations show that massive galax¬ 
ies form their stars early, and that the star formation in 
many of these massive objects is quenched early, so that 
their stars evolve largely passively. There is some tension 
in the ability of models to produce enough massive galax¬ 
ies at early times (z it 2), and a dearth of very rapidly 
star forming objects observed in the sub -mm and FIR 
(ISomerville et ai1l2012l : iNiemi et al.ll20lijf . However, the 
evolution of the number of massive “quenched” galaxies 
in models with AGN feedback seems to match observa¬ 
tions reasonably well dKimm et al] 120091 : iBrennan et al.l 
120151b 

Low mass galaxies seem to present a more thorny 
set of problems, which we refer to collectively as the 
“dwarf galaxy conundrum”. Models that reproduce the 
low-mass end of the observed stellar mass function lo¬ 
cally, generically overproduce low-mass (m* 10 10 Mq) 

galaxies at redshifts 0.5 it z it 2. Moreover, low-mass 
galaxies apparently have (specific) star formation rates 
that are too low over the same redshift range. The 
stellar ages predicted by our models for these galax¬ 
ies are too old compared to those derived for nearby 
galaxies based on ‘archaeological’ evidence. A summary 
of these problems, demonstrated for several indepen- 
dently developed sem i -analytic models, was p resented in 
iFontanot et alJ (l2009l b lWeinmann et al] (l2012l l presented 
a similar study that showed that the same problems also 
occur in numerical hydrodynam ic simulations, and re¬ 
cently ISomerville fe Pavel (|2014l ) showed that the prob¬ 
lem persists to varying degrees in most state-of-the-art 
SAMs and cosmological liydrodynamical simulations. It 
has been suggested that these problems might be due 
to inaccurate recipes for star formation, and that they 
might be cured by implementing metallicity dependent 
recipes for H 2 formation and H 2 -based star formation 
(iKrumholz fc Dekelll20 l j : iKuhlen et, alJl20l3 b This was 
one of the original motivations for the work we present 
here. 

The purpose of this paper is to present the details 
of how we incorporate partitioning of gas into an atomic, 
molecular and (optionally) ionized component in our ex¬ 
isting semi-analytic models, how we self-consistently im¬ 
plement an H 2 -based star formation recipe, and how sen¬ 
sitive our results are to details of the implementation. 
We explore three different recipes for the partitioning 
of gas into di fferent phases: the pressure-based recipe 
of lBIitz fe Rosolowskvl d2006l . BR) and two metallicity- 
b ased recipes, that of Krum holz et al. (KMT) and that 
of lGnedin fe Kravtsovl (l201ll . GK). We compare the pre¬ 
dictions of these three new models with those using the 
“classic” Kennicutt-Schmidt (KS) star formation recipe 
with no gas partitioning. In addition, we explore several 
different empirically motivated H 2 -based star formation 
recipes. 

This paper is pa rt of a series of related works. In 
IPopping et al.l (l2014d . PST14), we presented predictions 
for the atomic and molecular gas content of galaxies, 
and its evolution with redshi ft from z = 6- 0 , using the 
same models presented here. IPopping et alJ (l2014bh ex- 
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tended these models by carrying out radiative transfer 
calculations to predict sub-mm line emission luminosities 
from several atomic and molecu lar s pecies , including CO, 
HCN, C + , and [OI]. In lBerrv et all J20l4l . we presented 
predictions for the properties of objects that would be se¬ 
lected as Damped Lyman-a systems (DLAS) in absorp¬ 
tion against background quasars, again using the same 
model framework described here. In this paper, we focus 
on quantities pertaining to the stellar content, SFR, and 
metal content of galaxies and their evolution since z ~ 6. 
I 11 addition, we explore a wider variety of model variants 
than presented in the earlier works. 

The outline of the paper is as follows. In 92. II we out¬ 
line the basic framework of the semi-analytic models and 
the treatment of structure formation, gas cooling and in¬ 
fall, chemical evolution, and starbursts and morphologi¬ 
cal transformation via galaxy mergers. In 92.2l we describe 
our approaches for partitioning cold gas into an atomic, 
molecular, and (optionally) ionized component, in 92.3l we 
describe the new H 2 -based star formation recipes, and in 
Ea we describe our implementation of metal enhanced 
winds. In ED we describe how we choose the values of 
the free parameters in our models, and summarize their 
values. In 13. 1 1 we show how the star formation histories 
and build-up of stars, gas, and metals as a function of 
halo mass are impacted by the different recipes for gas 
partioning and star formation, and other details of our 
model implementation. In 93.21 we show predictions for 
the relationship between total gas density and SFR den¬ 
sity in our models. I 11 93.31 we present predictions for the 
stellar mass functions and stellar fractions, specific star 
formation rates, gas depletion timescales, and gas and 
stellar phase metallicities over cosmic time from z ~ 6 to 
the present. We discuss our results in m and summarize 
and conclude in El 


2 MODEL DESCRIPTION 

The semi-analytic models us ed here have been 
descri be d in detail in ISomerville &; Primackl 
(1 199(4). ISomerville et a 1 l200il 'l and most recently 
in ISom erville et all ( 2008al . hereafter S08) and 
ISomerville et all ( 20121 . S12). The Santa Cruz mod- 
eling framew ork has also recently been described in 
IPorter et all 12014) . We refer the reader to those papers 
for details. 


2.1 The Semi-Analytic Model Framework 

This section describes the aspects of the semi-analytic 
models that have been documented in previous papers. 
Therefore we give a relatively brief description of these 
ingredients here. 

In this work, the merging histories (or merger trees) 
of dark matter haloes are constructed based on the 
Extended Press-Sche chter (EP S) formalis m using the 
method described in ISomerville fc Kolattl d 19991 3 . with 
improvements described in S08. These merger trees 
record the growth of dark matter haloes via merging and 
accretion, with each “branch” representing a merger of 
two or more haloes. We follow each branch back in time 


to a minimum progenitor mass M ies , which we refer to 
as the mass resolution of our simulation. Our SAMs give 
nearly identical results when run on the EPS merger trees 
or on merger trees extracted from dissipationles s N-body 
simulations iLu et all l20l4 : IPorter et al.l l20l4) . We use 
EPS merger trees here because they allow us to attain ex¬ 
tremely high resolution, which is important for this study. 
We resolve halos down to M res = 10 10 M@ for all root ha¬ 
los, and below root halo masses of M res = 10 10 Mq, we 
set M res = 0.01 M roo t, where M roo t is the mass of the root 
halo. Our root halos cover a range from Mh = 5 x 10 8 Mq 
to 5 x 10 14 M 0 . 

When dark matter haloes merge, the central galaxy 
of the largest progenitor becomes the new central galaxy, 
and all others become ‘satellites’. Satellite galaxies lose 
angular momentum due to dynamical friction as they or¬ 
bit and may eventually merge with the central galaxy. 
To estimate this merger timescal e we use a variant of 
the C handrasekhar formula from iBovlan-Kolchin et al.l 
(|2008l l. Tidal stripping and destruction of satellites are 
also modeled as described in SOS. 

Before the Universe is reionised, each halo contains 
a mass of hot gas equal to the universal baryon fraction 
times the virial mass of the halo. After reionisation, the 
photo-ionising background suppresses the collapse of gas 
into low -mass haloes. We use the fitting functio ns pro¬ 
vided bv lOnedinl (12001)! ') and lKravtsov et all (l2004 ~). based 
on their hydrodynamic simulations, to model the fraction 
of baryons that can collapse into haloes of a given mass 
after reionisation, and assume that the universe was fully 
reionized by a = 11. 

When a dark matter halo collapses, or experiences 
a merger that at least doubles the mass of the largest 
progenitor, the hot gas is assumed to be shock-heated to 
the virial temperature of the new halo. This radiating 
gas then gradually cools and collapses. The cooling rate 
is estimated using a simple spherically sym metric model 
simila r to the one originally suggested bv I White fc Fren^ 
(Il99lh . Details are provided in S08. 

We assume here that the cold gas is accreted only by 
the central galaxy of the halo, although in reality satel¬ 
lite galaxies probably also continue to accrete some cold 
gas after they cross the virial radius of their host. I 11 
addition, we assume that all newly cooling gas initially 
collapses to form a rotationally supported disc. The scale 
radius of the disc is computed based on the initial angular 
momentum of the gas and the halo profile, assuming that 
angular momentum is conserved and that the self-gravity 
of the collapsing baryons causes contraction of the mat¬ 
ter in the inner pa r t of the halo dBlumenthal et all [l986 ; 
iFlores et al.l Il993l : iMo et ahl Il99sh . This approach has 
been shown to reproduce the observed size versus stellar 
m ass relation for d i sc-dom inated galaxies from z ~ 0- 
2 d Somerville et all 12008b '). In PST14 we also showed 
that our models reproduce the sizes of Hi disks in nearby 
galaxies, and sizes of CO disks out to z ~ 2. 

Star formation occurs in two modes, a normal “disc” 
mode in isolated discs, and a merger-driven “starburst” 
mode. Star formation in the disc mode is modelled as de¬ 
scribed in Section[233]below. The efficiency and timescale 
of the merger driven burst mode is modeled as described 
in S08 and is a function of the merger mass ratio and the 
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gas fractions of the progenitors. The treatment of merger- 
driven bursts is based on the results of hydrodynamic 
simu l ations of binary gala xy mergers dRobertson et all 
l200fJ : iHonkins et al.ll2009al f . 

Some of the energy from supernovae and massive 
stars is assumed to be deposited in the ISM, resulting 
in the driving of a large-scale outflow of cold gas from 
the galaxy. The mass outflow rate is 

( Vo \ arh 

rh 0 ut = £sn ( y- J rh * (2) 

where V c is the maximum circular velocity of the galaxy 
(here approximated by Vhiax of the dark matter halo), rh* 
is the star formation rate, csn and «sn are free param¬ 
eters, and Vo = 200 knr/s is an arbitrary normalization 
constant. Some fraction of this ejected gas escapes from 
the potential of the dark matter halo, while some is de¬ 
posited in the hot gas reservoir within the halo, where it 
becomes eligible to cool again. The fraction of gas that 
is ejected from the disc but retained in the halo, versus 
ejected from the disc and halo, is a function of the halo 
circular velocity (see S08 for details), such that low-mass 
haloes lose a larger fraction of their gas. 

The gas that is ejected from the halo is kept in a 
larger reservoir, along with the gas that has been pre¬ 
vented from falling in due to the photoionizing back¬ 
ground. This gas is assumed to accrete onto the halo on a 
timescale that is proportional to the halo dynamical time 
(see S08 for details). 

Each generation of stars produces heavy elements, 
and chemical enrichment is modelled in a simplified man¬ 
ner using the instantaneous recycling approximation. For 
each parcel of new stars dm*, we also create a mass of 
metals dMz = ydm*, which we assume to be instanta¬ 
neously mixed with the cold gas in the disc. The yield 
y is assumed to be constant, and is treated as a free pa¬ 
rameter. When gas is removed from the disc by supernova 
driven winds as described above, a corresponding propor¬ 
tion of metals is also removed and deposited either in the 
hot gas or outside the halo, following the same propor¬ 
tions as the ejected gas. Ejected metals also “re-accrete” 
into the halo along with the ejected gas, as described 
above. 

Mergers are assumed to remove angular momentum 
from the disc stars and to build up a spheriod. The 
efficiency of disc destruction and spheroid growth is a 
function of progenitor gas fraction and merger mass ra¬ 
tio, and is parameterized b ased o n hydrodynamic simu¬ 
lations of disc-disc mergers llHopkins et al.|[2009al i. These 
simulations indicate that more “major” (closer to equal 
mass ratio) and more gas-poor mergers are more efficient 
at removing angular momentum, destroying discs, and 
building spheroids. Note that the treatment of spheroid 
formation in mergers used here has been updated rela- 
tive to S08 as des cribed in iHopkins et al.l (I2009bh and 
IPorter et al.l (l2014h . We do not include a disk instability 
driven mode for spheroid growth in the models presented 
here. 

In addition, mergers drive gas into galactic nuclei, 
fueling black hole growth. Every galaxy is born with a 
small “seed” black hole (BH; here we adopt M see a ~ 
1.0 x 10 4 Mq). Following a merger, any pre-existing black 


holes are assumed to merge immediately, and the result¬ 
ing hole grows at its Eddington rate until the energy be¬ 
ing deposited into the ISM in the central region of the 
galaxy is sufficient to significantly offset and eventually 
halt accretion via a pressure-driven outflow. This results 
in self-regulated accretion that leaves behind black holes 
that naturally obey the observed correlation between BH 
mass and spheroid mass or velocity dispersion. Our mod¬ 
els produce good agreement with the observed luminosity 
function of X-ray /optical/ IR detected quasars and AGN 
dHirschmann et al.l 1 20 12bl V . 

A second mode of black hole growth, termed “radio 
mode”, is associated with powerful jets observed at ra¬ 
dio frequencies. In contrast to the merger-triggered, ra- 
diatively efficient mode of BH growth described above 
(sometimes called “bright mode” or “quasar mode”), in 
which the BH accretion is fueled by cold gas in the nu¬ 
cleus, here, hot halo gas is assum ed to be accr eted accord¬ 
ing to the Bondi-Hoyle model |Bon dl l 195211 . This leads 
to accretion rates that are typically < 10“ 3 times the Ed¬ 
dington rate, so that most of the BH’s mass is acquired 
during episodes of “bright mode” accretion. However, the 
radio jets are assumed to couple very efficiently with the 
hot halo gas, and to provide a heating term that can par¬ 
tially or completely offset cooling during the “hot flow” 
mode (we assume that the jets cannot couple efficiently 
to the cold, dense gas in the infall-limited or cold flow 
regime). 

2.2 Multi-phase Gas Partitioning 

In this section, we describe in detail the updates to the 
model ingredients that are explored in this paper. These 
include partitioning of the cold gas in galactic disks into 
an ionized (Hu), atomic (Hi), and molecular (H 2 ) com¬ 
ponent, an option to include metal enhanced winds, and 
a set of new H 2 -based star formation recipes. 

At each timestep, we compute the scale radius of 
the cold gas disc using the angular momentum based ap¬ 
proach described above, and assume that the total (Hi 
+H 2 +H 11 ) radial cold gas distribution is described by 
an exponential with scale radius r gaa . We do not attempt 
to track the scale radius of the stellar disk separately, but 
make the simple assumption that r, = x gas r gas , with Xgas 
fixed to match obser ved stellar scale lengths at 2 = 0. 
iBigiel fc Blitzl (120121 ') showed that this is a fairly good 
represen tation, on avera ge, of the disks of nearby spirals 
(see also lKravtso\l20l3 l. We then divide the gas disk into 
radial annuli and compute the fraction of molecular gas, 
/h2(i~) = SH 2 (r)/[SH 2 (r) + Em(r)], in each annulus, as 
described below. We use a fifth order Runga-Kutta inte¬ 
gration scheme to compute the integrated mass of Hi and 
H 2 in the disk, and the integrated SFR, at each timestep. 

2.2.1 Ionized gas associated with galaxies 

Most previous semi-analytic models have neglected the 
ionized gas associated with galaxies, which may be ion¬ 
ized either by an external background or by the radiation 
field from stars within the galaxy. Here we inves tigate a 
simple analytic model motivated by the work of IGnedirJ 
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J2012I) . We assume that some fraction of the total cold 
gas in the galaxy, /ion,int, is ionized by the galaxy’s own 
stars. In addition, a slab of gas on each side of the disk 
is ionized by the external background radiation field. As¬ 
suming that all gas with a surface density below some 
critical value Ehii is ionized, we have 


Shii 

^T 


1 + In 



+ 0.5 In 



2 ' 


where So = m c oid/(27rr gaa ) 2 is the central surface den¬ 
sity of the cold gas (m co id is the mass of all cold gas in 
the disk and r gaa is the scale radius of the gas disk). We 
typically assume /i on ,mt = 0.2 (as in the Milky Way) and 
Ehii = 0.4Mgpc -2 , as in IGnedinl d2012l '). Applying this 
model within our SAM gives remarkably good agreement 
with the ionized fraction s as a f u nction of circular ve¬ 
locity shown in Fig. 2 of IGnedinl (|2012l l. obtained from 
hydrodynamic simulations with time dependent and spa¬ 
tially variable 3D radiative transfer of ionizing radiation 
from local sources and the cosmic background. 


2.2.2 Molecular gas: pressure based partitioning 


We consider three approaches for computing the molec¬ 
ular gas fractions in galaxies. The first is based 
on the empirical pressur e-based recipe presented by 
iBlitz fe Rosolowskvi (20 06jf. and will b e referred to as the 
BR recipe. Blitz fc Rosolowskvi d2006l 'l found a power-law 
relation between the disc mid-plane pressure and the ra¬ 
tio between molecular and atomic hydrogen, i.e., 





(3) 


where Eh 2 and Shi are the H 2 and Hi surface den¬ 
sity, Po and qbr are free parameters that are obtained 
from a fit to the observational data, and P m is the mid¬ 
plane pressure acting on the galactic disc. We adopted 
logPo/pB = 4.23 cm 3 K and o br = 0.8 based on obser¬ 
vations from iLerov et all d2008lf . 

We estimate the hydrostatic pressure as a function 
of the distance from the center of the disk r as 

P(r) = |GS gaa (r)[E gaa (r) + / a (r)S.(r)] (4) 

where G is the gravitational constant, S gas is the cold gas 
surface density, E* is the stellar surface density, and f a 
is the ratio of the vertical velocity dispersion s of th e gas 
and stars: / CT (r) = ^ 21 . Following iFu et al. I (1201011 , we 

adopt /o-(r) = O.l-y/St.o/S,, where E„,o = m*/{2rrl), 
based on empirical scalings for nearby disk galaxies. 


2.2.3 Molecular gas: metallicity based partitioning 

Gnedin and Kravtsov (jOnedin fc Kravtsov! i201di . 2011 ! 
performed high-resolution “zoom-in” cosmological simu- 
la tions with the Adapti ve Refinement Tree (ART) code 
of lKravtsov et all (Il997l l , including gravity, hydrodynam¬ 
ics, non-equilibrium chemistry, and simplified on-the-fly 
radiative transfer. These simulations are therefore able to 
follow the formation of molecular hydrogen through pri¬ 
mordial channels and on dust grains, as well as dissocia¬ 
tion of molecular hydrogen and self- and dust- shielding. 


1 Gnedin fe Kravtsovl (l201ll f presented a fitting func¬ 
tion based on their simulations, which effectively param¬ 
eterizes the fraction of molecular hydrogen as a function 
of the dust-to-gas ratio relative to the Milky Way, Dmw, 
the UV ionizing background relative to the Milky Way 
Fmwi and the neutral gas surface density S hi+h 2 - The 
fraction of molecular hydrogen is given by 


fn 2 


1 + 


S» 

S hi+h 2 


-2 


where 
S* = 
A = 
9 = 

s = 

a — 
D * = 


2 A 4 / 7 1 

20Mq P c -2 —- 

D MW yj 1 —(— Umw D 

ln(l + 5 P^(C7 M w/15) 4/7 ) 


2 

MW 


1 + as + s 2 
1 + s 
0.04 


D * + Umw 


_ Pmw/2 
1 + (Pmw/2) 2 

1.5 x 10 -3 ln(l + (3Pmw) 17 ) 


We take the dust-to-gas ratio to be proportional to 
the cold gas phase metallicity in solar units Bmw = 
Z/Zq. The local UV background relative to the MW 
is assumed to scale in proportion with the global 
SFR of the galaxy in the previous time step rela¬ 
tive to the MW SFR, Umw = , where we 

c hoose SFRmw = 1.0 M g yr - 1 (IMurrav fe Rahmanl 

I 2 OI 0 I ; [Robitaille fe Whitnevl I 2 OI 0 II . We refer to this as 
our ‘fiducial’ GK model. We also investigate the results 
of keeping Umw fixed to the Milky Way value, which we 
refer to as the GKFUV model (fixed UV field). 

An alternate approach based on similar physi¬ 
cal processes was presented in a series of papers 
by Kr um holz and collaborators llKrumholz et al.l 120081 . 
l2009allbh . iKrumholz et al.l ll2009blf developed an analytic 
model for the molecular fraction in galaxies, based on 
the ansatz that the interplay between the interstellar ra¬ 
diation field and molecular self-shielding determines the 
molecular fraction. They presented a fitting function: 


}h 2 =1- 


1 + 


3 s \ 
41 + Sj 



where s = ln(l + 0.6x)/(0.04E comp ,oZ'), X = 0.77(1 + 
3.1Z' 0 - 365 ), S = 0.0712(0.Is -1 + 0.675) -2 ' 8 , S comp , 0 = 
S comp /(lM 0 pc -2 ), and Z' = Z/Zq. The quantity de¬ 
noted Ecom P is the surface density wit hin a ~ 100 pc 
sized a tomic-molecular cloud complex. IKrumholz et al.l 
(|2009bh suggest using a “clumping factor” to apply the 
model to simulations with spatial resolution coarser than 
100 pc; i.e. S com p cE H i+h 2 , where c > 1, and E H i+h 2 
is the neutral gas surface density on some larger scale. 
The appropriate value of c depends on this spatial scale, 
where c —> 1 as the scale approaches 100 pc. We refer to 

this as the KMT gas partitioning r ecipe._ 

Both _ iGnedin &; Kravtsovl (1201 ll ! and 

IKrumholz et al. ( 2009bi l note that the fitting func¬ 
tions, as well as perhaps (in the case of KMT) the 
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underlying assumptions of the model, begin to break 
down at metallicities lower than about l/50th of the 
Solar value. 

The KMT and GK fitting functions above character¬ 
ize the formation of H 2 on dust grains, which is the dom¬ 
inant mechanism once the gas is enriched to more than 
a few hundredths of Solar metallicity. Other channels for 
the formation of H 2 in primordial gas must be responsible 
for producing the molecular hydrogen out of which the 
first stars were formed. Hydrodynamic simulations con¬ 
taining detailed chemical networks and analytic calcula¬ 
tions have shown that H 2 can form in metal-free gas in 
dark matter halos ab ove a critical mass M CI jt ~ 10 5 M@ 
(e.g., iNakamura fc Umemurai i200ll ; [Gloveil l2013j l. This 
gas can then form “Pop III” stars w hich can enrich 
the surrounding ISM to Zm ~ 10~ 3 Z 0 (I Schneider et al.l 
120021 : [Greif et all20ld : IWise et al.ll2012l l. These processes 
take place in halos much smaller than our resolution limit. 
We represent them by setting a “floor” to the molecular 
hydrogen fraction in our halos, /h 2 ,floor- In addition, we 
“pre-enrich” the initial hot gas in halos, and the gas ac¬ 
creted onto halos due to cosmological infall, to a metal¬ 
licity of Z pre _ enric h. We adopt typic al values of / H 2 .floor = 
10~ 4 and Z pre —enrich = 10 -3 Z 0 (1 11 aim an et al. 1 19961 : 
iBromm &; Larsonl l2004l l. We explore the sensitivity of 
our results to these parameters in Section 13.11 Note that 
observations of resolved stars in the halo of our Galaxy 
and local dwarfs hav e revealed stars with metallicities 
below Z ~ 10 -3 Zm llTolstov et ahl 120091 : ISchorck et al.l 
120091 : iKirbv et al.l 1201 lh . precluding much higher values 
for Zp re _ enr i c h. 

2.3 Star Formation Recipes 

2.3.1 The Kennicutt-Schmidt (KS) Recipe 

The KS recipe (K ennicu iili'H i assumes that the surface 
density of star formation in a galaxy is a function of the 
total surface density of the cold neutral gas (atomic and 
molecular), above some threshold surface density E cr it. 

The star formation rate density (per unit area) for 
Egas > E C rit is given by: 

Esfr = Asf Eg as JVsF , (5) 

where Esfr = 0 for E gaa < E cr it- This recipe is the same 
one used in most of our previously published SAMs (S08, 
S12), and is similar to recipes commonly adopted in many 
other SAMs. The values of the free parameters are given 
in Table |T] 


2.3.2 Molecular Hydrogen-based Recipes 

In the same spirit as the KS recipe, we use empirical 
relation ships from observat ions to motivate our H 2 -based 
recipes. iBigiel et al.l (|2008T) found, based on observations 
of spiral galaxies from the THINGS survey, that the star- 
formation rate surface density can be directly related to 
the surface density of molecular gas, i.e. 

Esfr= (iOMopc-0 (6) 



Figure 1. Empirically based star formation recipes used as 
input in our models, and from the literature. The solid purple 
line shows the two-slope H 2 -based recipe (Big2); the dashed 
lavender line shows the single slope H 2 -based recipe (Bigl); the 
dot-da shed dark red line shows the recipe based on the anal¬ 
ysis of I Narayanan et al.1 (|2012h : the dotted green line shows 
the (Hi + H 2 ) based “classic” KS recipe. The vertical dashed 
line shows the critical total gas surface density used in our 
models that implement the classic KS recipe (£ C rit )- The star 
symbol shows the relation derived by lSharo n et al l d2013h for 
an ext reme starburst galaxy. Note that the [Naravana net al.l 
d2012T) results are shown for reference only; we do not show 
the results of incorporating this recipe in our models here. 

with Nsf — 1 (see also IBigiel et al.l l201ll ; I Leroy et al.l 
l2013lb Observations of higher density environments sug¬ 
gest that above some critical H 2 surface density, the 
slope of the relation described in Eqn. |S] steepens 
dNaravanan et ahl 120121 '). We therefore also consider a 
two-part scaling law given by: 

Esfr = ^ SF (1OM0PC -2 ) ( x + sfi) 

The values of the parameters Asf, Nsf, and EH 2 ,crit are 
given in Table |T] 

A star formation relation that changes slope above 
a critical density is also expecte d based on theoretical 
grounds. iKrumholz et al.l (I2009bl) adopt the star forma¬ 
tion relation: 

/ y \ n sf 

Esfr = AsfEr 2 ( gas ) (8) 

where Nsf = —0.33 for Eg a a/E C rit < 1 and Nsf = 0.33 
for Egas/Ecrit > 1. KMT adopt Asf = 1/2.6 Gyr -1 and 
E cr it = 85 M 0 pc~ 2 . We adopt the same parameter values 
in our “KMT” model. 

A comparison of various star formation relations in 
the literature, and used in this work, is shown in Fig. [T] 

2.4 Metal Enhanced Winds 

Most of our previous models have assumed that metals 
are ejected from galaxies with the same efficiency as the 
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gas, i.e. with the same mass loading factor r/ = rhout/m*. 
However, since metals are produced by the same massive 
stars and supernovae that are believed to drive galac¬ 
tic outflows, it is possible that metals are preferentially 
ejected (i.e., have a higher effective mass loading fac¬ 
tor than the gas averaged over the whole disk). Because 
two of our recipes for gas partitioning depend on the gas 
metallicity, the dispersal of metals in our models has a 
potentially important impact on our results. We therefore 
include an optional treatment of metal-enhanced winds 
in our models. 

We base our parameter ization of metal-enhance d 
winds on the approach used in lKrumholz & Dekell ( 20121 ). 
in part because we want to be able to compare our re¬ 
sults with theirs. The fraction of metals that is ejected is 
parameterized by: 

C = Cl°exp (-Mh/M r et) 

where both (ji 0 and M ret are free parameters, and Mh is 
the virial mass of the halo. The modified equation for the 
evolution of the mass in metals in the cold gas phase is 
then: 

— y( 1 R)(l ^)t? 2* -f ZUiotUZinf ^coldUlout 

where R is the recycled fraction, y is the chemical yield, 
Zhot is the metallicity of the hot gas, Z co id is the metal¬ 
licity of the cold gas, and m„, mint, and m ou t are the 
star formation rate, inflow rate of gas from the hot halo 
into the disk, and the outflow rate of gas from the disk, 
respectively. 


2.5 Calibrating the Free Parameters 


As in any cosmological simulation, we must parameter¬ 
ize the sub-grid physics in our models. In keeping with 
common practice, we choose the values of the free param¬ 
eters by tuning to a subset of observations in the local 
universe. In this subsection we summarize the values of 
the free parameters used here (see Table []J and the ob¬ 
servations we used to constrain them. The parameters 
are t h e same as t hos e used in the models presented in 
IPopping et al.l (l2014d ~). 

We assume values for the cosmological parameters 
consistent with the five year WMAP results (WMAP5): 


fl m = 0.28, Pa = 0.72, H o = 70.0, as = 0.81, and n s = 
0.96 l|Komatsu et al.ff2009l ). We note that these values are 
generally consistent with those obtained from the analy- 
sis of the seven-year WMAP data release llKomatsu et al.l 
120101 '). The adopted baryon fraction is 0.1658. We assume 
a recycled fraction of R = 0.43, as approp riate for a 
Chabrier stellar Initial Mass Function jchnbrier} 20031). 

As discussed in S08 (see also IWhite et al. 2014 ), in 
our models the supernova feedback parameters mainly 
control the low-mass end of the stellar mass function 
(m* M c har, where M c h ar is the characteristic mass of 
the “knee” in the Schechter function describing the stel¬ 
lar mass function), or equivalently, the fraction of baryons 
that is turned into stars in halos with Mh S 10 12 Mq. On 
the other side, the efficiency of the “radio mode” AGN 
feedback (one can think of this schematically as the effi¬ 
ciency with which radio jets couple to and heat the hot 


intragroup and intracluster medium) controls the number 
density of massive galaxies m» M c h ar , or the fraction 
of baryons that are able to turn into stars in massive 
halos (Mh it 10 12 Mq). We tune the parameters control¬ 
ling supernova feedback and AGN feedback to reproduce 
the observed stellar mass function at z = 0 in our tradi¬ 
tional KS model (see SOS for details). These parameters 
are then kept fixed as we explore the effects of varying 
the modeling of gas partitioning and star formation. 

The parameters of the star formation recipe mainly 
control the fraction of cold gas in galaxies, and do 
not strongly affec t the z = 0 stellar mass function 
(IWhite et al.l I2014| j . We require the parameters charac¬ 
terizing our star formation recipes to lie within the obser¬ 
vational uncertainties from recent empirical constraints, 
and tune them within these limits to match the total gas 
fractions as a function of galaxy stellar mass at z = 0 
(see PST14). 

The chemical yield y could in principle be obtained 
from stellar evolution models, but these model yields are 
uncertain by a factor of ~ 2, and the single-element 
instantaneous recycling approach to chemical evolution 
that we are using here is somewhat crude, so we instead 
treat the chemical yield as a free parameter (though we 
restrict it to be in the expected range). We tune our yield 
to match the normalizati on of the observed stel lar metal¬ 
licity vs. mass relation of lGallazzi et al.l (|2005l ). 

We take the parameter values for the BR gas 
partitioning recip e from the observational results of 
iLerov et al! |2008|). We impl e ment the GK recipe as it is 
given in Gnedin fc Kravtsov! J201 ill , with no tunable pa¬ 
rameters. The KMT recipe has one free parameter, the 
clumping factor of the g as. We adopt c = 5, following 
iKrumholz &; Dekell (12012 ). 

Our new models predict the fraction of cold gas 
in different phases: ionized, atomic, and molecular. We 
showed the predictions of our two fiducial models (GK 
and BR) for the fraction of Hi and H 2 as a function of 
galaxy internal stellar density and stellar mass in PST14. 
It is encouraging that our new models reproduce these ob¬ 
served scalings for nearby galaxies quite well without any 
additional tuning. We emphasize that the only new free 
parameters in the gas partioning recipe have been taken 
directly from observations (in the case of the BR recipe) 
or from numerical simulations (in the case of KMT and 
GK). 


3 RESULTS 

3.1 Effects of Varying Model Ingredients, 
Parameter Values, and Resolution 

In this subsection we explore the sensitivity of our model 
results to our new model ingredients and parameter val¬ 
ues related to gas partitioning and star formation, as 
well as to our numerical resolution. To illustrate these 
effects, we show the properties of the largest progenitor 
galaxy as a function of time (or redshift) in a set of halos 
with masses at z = 0 ranging from logM^/Mo = 10.0- 
11.5 (except in one case, where we show a more massive 
halo with log Mh/Mq = 12.5). The variations that we 
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Table 1. Summary of Model Parameters 


parameter 

description 

section defined 

value 

e SN 

supernova feedback efficiency 

0 

1.5 

0!SN 

supernova feedback slope 

m 

-2.2 

y 

chemical yield 

0 

1.6 Z 0 

« AGN 

radio mode AGN feedback 

S08 §2.11, Eqn. 20 

3.8 x 10 -3 

Xgas 

ratio of stellar to gas scale length 

m\ 

0.59 

Shii 

critical density for ionized gas 

12.2.11 

0.4 M 0 pc -2 

/ion,int 

internal fraction of ionized gas 

12 XT 1 

0.2 

Po 

pressure scaling in BR recipe 

12.2.21 

4.23 kg cm 3 K 

a BR 

slope in BR recipe 

12.2.21 

0.8 

C 

clumping factor in KMT recipe 

|2.2.3| 

5 

fl!2, floor 

primordial H 2 fraction 

12.2.31 

10" 4 

^pre—enrich 

metallicity due to Pop III stars 

|2.2.3| 

10“ 3 

Clo 

metal enhanced winds normalization 

EU 

0.1 


metal enhanced winds mass scale 

I2~T1 

0.9 

A sf (KS) 

star formation relation normalization 

12.31 Eqn.0 

1.1 x 10 -4 

Nsf (KS) 

star formation relation slope 

12731 Eqn. 0 

1.4 

Scrit (KS) 

critical density for SF 

12731 Eqn. 0 

6 M 0 pc~ 2 

Asf (Bigl, Big2) star formation relation normalization 

12.31 Eqn. 0 [7] 

4.0 x 10“ 3 

TVsf (Bigl, Big2) star formation relation slope 

1231 Eqn. 0 0 

1.0 

^H2,crit 

critical H 2 density 

12731 Eqn. [7] 

70 M 0 pc~ 2 


Table 2. Summary of Model Variants 


model 

H 1 /H 2 partitioning 

SF law 

metal-enhanced winds 

KS “fiducial” 

none 

KS 

N 

BR “fiducial” 

BR 

Big2 

N 

GK “fiducial” 

GK, U M w °c SFR 

Big2 

N 

GK+Bigl 

GK 

Bigl 

N 

GKFUV 

GK, U M w = 1 

Big2 

N 

BR+Bigl 

BR 

Bigl 

N 

KMT+Bigl 

KMT 

Bigl 

N 

KMT 

KMT 

KMT 

N 

KMT+MEW 

KMT 

KMT 

Y 


explore here mainly affect lower mass halos, and produce 
no significant differences for galaxies in halos more mas¬ 
sive than those shown. We use the same merger trees for 
each model, and fix all model properties that are chosen 
from random distributions to their average values. Each 
panel shows the average over 60 different realizations of 
halos with the specified final mass. For each experiment 
we show the stellar mass, total neutral cold gas mass 
(Hi +H 2 ), SFR, H 2 fraction, and gas phase metallicity. 
To facilitate comparison, the stellar mass, gas mass, and 
SFR are normalized by dividing by the root halo mass at 
z = 0. 

As a first basic check, we test the impact of varying 
the mass resolution of our merger trees (Fig. 0. Note 
that what we mean by the ‘mass resolution’ here is the 
mass of the smallest halos that are tracked in the merger 
tree. This is not equivalent to the particle mass in an 
Wbody simulation, but rather to the smallest halo mass 
that can be robustly identified. We see from this test that 


in order to robustly reconstruct the whole halo mass ac¬ 
cretion history back to z ~ 10, we require a minimum 
halo mass of ~ 1/100 of the root mass at the output red- 
shift. Accordingly, we impose this condition on all halos 
in the runs used in this work. It is reassuring to see that, 
once the halo mass accretion history is well resolved, our 
SAM predictions converge extremely well (note that we 
do not retune the free parameters when we change the 
mass resolution). 

In Fig. [3] we test for possible sensitivity to the val¬ 
ues of two parameters that we introduce to simulate the 
formation of stars in primordial gas, the metallicity of 
the “pre-enriched” gas Z pre - e nr i c h (most relevant for the 
metallicity-based recipes), and the primordial molecu¬ 
lar hydrogen fraction /h 2 ,floor (most relevant for the BR 
recipe). Leaving all other settings of the fiducial GK 
model fixed, we vary Zpre-enrich from its fiducial value 
of ICG 3 by one order of magnitude downwards, to 10 -4 , 
and upwards to 0.01. Overall, the impact of even such 
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redshift redshift redshift redshift 



r * , , / , , ;, /. , • . 

0.5 1 2 4 6 10 0.5 1 2 4610 0.5 1 2 4610 0.5 1 2 4610 

time (Gyr) time (Gyr) time (Gyr) time (Gyr) 


Figure 2. From top to bottom, colored lines show the stellar mass, cold neutral gas mass (Hi +H2), and SFR normalized by 
the mass of the root halo at z = 0 for the largest progenitor galaxy as a function of cosmic time (redshift). The H 2 fraction 
(/h 2 = rr> H 2 /( rn H2 + rn Hl ) and gas phase metallicity (in solar units) are also shown. Gray lines show the maximum baryon 
fraction in the halo, /), A/), (1), where M^(t) is the mass of the largest progenitor halo at time t and /(, is the universal baryon 
fraction (different gray lines correspond to different resolutions; lower resolution runs cannot resolve the halo mass accretion history 
as far back in time). In this experiment, we test the dependence of our results on the mass resolution of our merger trees, varying 
the mass resolution by two orders of magnitude. Results are shown for a mass resolution of 10 10 , Mg (solid purple), 10 9 Mg 
(dashed orange), and 10 s Mg (dotted green). The mass accretion histories are well-resolved when the mass resolution is at least 
1/100 the mass of the root halo. 


extreme variations is fairly minor. The most noticable 
impact is on the stellar and gas phase metallicities. The 
metallicity builds up earlier in models with higher val¬ 
ues of Zpre-enrich, as expected. As a result, the H 2 frac¬ 
tion is higher at earlier times in the model with higher 
Zp re _eiirich, leading to higher star formation efficiency 
(SFE) and slightly lower gas fractions. Similarly, we var¬ 
ied the value of /h 2,floor from its fiducial value of 10 -4 
up and down by an order of magnitude. This has no dis- 
cernable effect on our results, and we therefore omit the 
corresponding figure. 


In a related experiment, we run our (otherwise) fidu¬ 
cial GK model with metal-enriched winds, described in 
Section GO Here, the metallicity of stellar driven out¬ 
flows can be metal-enhanced relative to the ISM by 
a factor that depends on the halo mass. As seen in 
Fig. H we find that metal enhanced winds can signif¬ 
icantly delay the formation of H 2 and stars in very 
low-mass halos (log Mh/Mq = 10.0), and cause the 
build-up of slightly more cold gas in low-mass halos 
(logAVMg ^ 10.5); note however that these halos host 
galaxies that are well below the detection limits of most 
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redshift redshift redshift redshift 



0.5 1 2 4 6 10 0.5 1 2 4610 0.5 1 2 4610 0.5 1 2 4610 
time (Gyr) time (Gyr) time (Gyr) time (Gyr) 


Figure 3. Same as Fig. [2] except here we compare the results of our fiducial GK model with different values for the “pre-enriched” 
gas metallicity (Z pre _ enr ; c h)i as shown in the key. Our results are quite insensitive to the value of this parameter within a reasonable 
range. 


surveys except in the very nearby Universe (m* ~ 10 7 - 
10 s Mg). Metal-enriched winds can also delay the build¬ 
up of metal-enriched gas even in more massive halos 
(log M h /M @ < 11.5). 

A new ingredient we have introduced into our mod¬ 
els is the tracking of gas that is photoionized either by 
an external radiation field or by internal sources. This 
gas is not eligible to form H 2 or stars. In Fig. [5] we show 
the galaxy properties in the fiducial GK model with and 
without tracking of Hu. Although our model predicts that 
ga laxies contain a subs tantial amount of Hu (see Fig. 2 
in [Popping et al]|2014df . partitioning this gas into a sep¬ 
arate reservoir has a very weak effect on our results. The 
only noticable effect is slightly lower H 2 fractions at high 
redshift, particularly in the lowest mass halos. 

In the next experiment, shown in Fig. [6] we investi¬ 


gate several different recipes for converting cold molecu¬ 
lar gas into stars within the fiducial GK model. We con¬ 
sider two vari ants of t he empirical recipe based on the ob¬ 
servations bv lBigiel et al.l (120081 ). In addition, we consider 
the recipe proposed by KMT based on theoretical argu¬ 
ments (see Section [2.31 for details). Bigl refers to Eqn. [fj] 
and Big2 to Eqn. [7] In contrast to most of our other ex¬ 
periments, these variations have almost no discernable ef¬ 
fect on the low mass halos. However, the Big2 recipe leads 
to significantly earlier build-up of stellar mass, more ef¬ 
ficient star formation at high redshift, and earlier metal 
enrichment in massive halos (log Mh/Mq (i, 12.0). This 
is owing to the non-linear dependence of the star for¬ 
mation efficiency on H 2 density at high densities in this 
model, and the positive correlation between halo mass 
and galaxy surface density in our models. 
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Figure 4. Same as Fig. [2] except here we compare the results of our fiducial GK model with and without metal enhanced winds. 
Metal enhanced winds can delay enrichment and, to a lesser extent, H 2 and star formation in low mass halos. 


Next we experiment with changing the recipe for 
partitioning gas into H 2 (Fig. Q. All other ingredients 
are the same as our fiducial GK models. We show the 
pressure-based BR model as well as two alternate metal- 
licity based models. Recall that in the fiducial GK model, 
/h 2 depends on the total gas density and metallicity as 
well as the local UV radiation field (which we scale with 
the global galaxy SFR). In the GKFUV model, we re¬ 
move the UV radiation field dependence by using the 
Milky Way value in all galaxies. In the KMT model, /h 2 
depends only on total gas density and metallicity, and has 
different dependencies on these quantities than the GK 
model. The results of this experiment are quite interest¬ 
ing. The predictions of the BR and fiducial GK models 
are quite similar, although the GK model tends to predict 
higher gas masses, lower H 2 fractions, and lower metallic- 
ities at early times in the two lowest mass halo bins. The 


KMT and GKFUV models also pr oduce similar results, 
as exp ected based on the findings of lKrumholz fe Gnedinl 
(I2011B . who also showed the two models to be very simi¬ 
lar. However, the KMT and GKFUV models predict sig¬ 
nificantly suppressed H 2 formation, leading to lower star 
formation rates and stellar masses, reduced chemical en¬ 
richment, and higher gas fractions in the two lowest halo 
mass bins. This is because galaxies in low-mass halos tend 
to have lower metallicities but also lower SFR. Therefore 
in our fiducial GK models, the lower metallicity, which 
makes H 2 formation less efficient, is mitigated by the 
lower SFR, which leads to weaker photo-dissociation and 
relatively higher H 2 fractions. 

In our penultimate experiment (Fig. |8}, we compare 
our two new fiducial multiphase gas recipes with the 
“classic” Kennicutt-Schmidt recipe (see Section 12.311 , in 
which all cold gas above a fixed surface density is eligible 
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Figure 5. Same as Fig. [2] except here we compare the GK model with and without the partioning of ionized gas (Ho) into a 
separate reservoir. Our results are nearly unchanged whether or not we include the ionized gas component. 


for star formation. This KS reci pe has been used in many 
previous SAMs (e.g. S08, S12, [Porter et al.ll20f3l . The 
build-up of stellar mass is almost identical in all three 
models. However, the cold gas mass is highest in the GK 
model and tends to be lowest in the KS model. The H 2 
fraction is also lower at early times in the lowest halo 
mass bin in the GK model. Note that the H 2 fraction 
shown for the KS model has been computed using the 
BR recipe in post-processing, but this has no impact on 
the star formation in the model. Overall, the degree of 
similarity between the three model results is quite sur¬ 
prising, given the rather different physical premises on 
which they are based. We discuss possible reasons for 
this in Sj4] 

In our final experiment, shown in Fig. [9] we compare 
the evolution in our fiducial GK model with variants that 
include combinations of recipes that are similar to those 


used in several published models from the literature. For 
example, the BR+Bigl model treats gas partitioning and 
conversio n of H 2 to stars usin g similar recipes to the BR 
model of lLagos et ali d2011bll and the “Bigiel + H 2 pre¬ 
scription 2” of Fu et alT i 2012 1. The KMT+Bigl contains 
simil ar ingredients t o the “Krumholz + H 2 prescription 
1” of lFu et al.1 ()2012h . In the KMT+MEW model, we use 
the KMT recipes for both gas partioning and star for¬ 
mation, as well a s including metal enhance d winds, as 
in the models of lKrumholz_fc Dekell (I2012T ). Note that 
the KMT model of Lagos et al.l ( 201 lbh adopts the KMT 
recipes for both gas partioning and star formation, but 
does not adopt metal-enhanced winds, so does not cor¬ 
respond exactly to any of the cases shown here. How¬ 
ever, adopting these choices in our models yields results 
very similar to the KMT+Bigl model shown. We empha¬ 
size that many other aspects of our models differ from 
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Figure 6. Same as Fig. [2] except here we compare different recipes for converting molecular gas into stars within our fiducial 
GK models: Big2, Bigl, and KMT. Note the different range of halo masses from the other plots. Here, the strongest effect seen is 
the more efficient production of stars and earlier enrichment in massive halos in the Big2 model. This is owing to the non-linear 
dependence of the star formation efficiency on H 2 density at high densities in this model. 


those used by other SAMs in the literature, so these 
may not correspond to the actual predictions of those 
models. This exercise is intended to shed some light on 
the effect of choosing different recipes for gas partition¬ 
ing and conversion of H 2 into stars in a controlled envi¬ 
ronment where all other aspects of the models are held 
fixed. Most other SAMs to date that have attempted to 
track multi-phase gas with a metallicity-based approach 
have done so using the KMT recipe. Our experiment 
shows that this may result in more delayed star forma¬ 
tion and enrichmen t in low-mass halo s than our fidu¬ 
cial models predict. iKrumholz fe Dekell (l2012li addition¬ 
ally adopt strongly halo mass dependent metal-enhanced 
winds. These two effects together lead to strong suppres¬ 


sion of star formation in low-mass galaxies, particularly 
at early times. 


3.2 Star Formation Relations 

Fig. [10] shows the relationship between the total neu¬ 
tral cold gas surface density T,hi+h 2 an d star formation 
rate density Esfr in our two new fiducial models with 
multiphase gas partitioning. In our previous generation 
of models, galaxies had a deterministic relation between 
T,hi+h 2 and Esfr given by the assumed KS relation (as 
plotted in Fig. [TJ, and Esfr was set to zero below the 
critical gas surface density E cr it (also shown in Fig.Q]). In 
our new models, neutral gas is ‘partitioned’ into Hi and 
H 2 , and only H 2 is allowed to participate in star forma- 
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Figure 7. Same as Fig. [2] except here we compare our fiducial GK and BR models with the KMT recipe for gas partioning 
and the GK recipe for gas partitioning with a fixed UV background field. The KMT recipe and the GKFUV recipes, which both 
neglect the dependence of H 2 formation on the local UV radiation field, both predict much lower H 2 fractions in low mass halos, 
especially at high redshift, leading to later stellar mass assembly, slightly higher overall gas fractions, lower SFE, and later metal 
enrichment. 


tion. Therefore the value of Esfr at a given T,hi+h 2 has 
a “second parameter” dependence. This second parame¬ 
ter is metallicity in the case of the GK (and KMT, not 
shown) recipes and disk mid-plane pressure (stellar sur¬ 
face density E*, to first order) in the BR recipe. Fig. [10] 
shows the the average metallicity of the cold gas in each 
galaxy, where each dot shows one annulus with radius 
500 pc. In the GK model, galaxies with higher metallic¬ 
ity have a higher Esfr for a given E hi+h 2 , as expected. 
However, we also see a similar dependence on metallicity 
in the BR model, although in this case it is not directly 
input into the model. The reason for this apparent de¬ 
pendence on metallicity is simply that E* and gas phase 
metallicity are highly correlated. The curvature in Esfr- 


F,hi+h 2 at low gas surface densities in both models is in 
goo d agreement with observations of nea rby spiral galax¬ 
ies llLerov et al.ll2008l : iBigiel et a.l .11 20081) . 

3.3 Evolution of Galaxy Populations 

3.3.1 Stellar Mass Functions and Stellar Fractions 

In this sub-section we examine the evolution of galaxy 
populations, which are directly comparable with obser¬ 
vations. We consider three main model variants: the fidu¬ 
cial versions of the GK, BR and KS models (see Table[2]|. 
The KS model is the star formation reci pe used in previ¬ 
ously published Santa Cruz SA Ms (e.g. ISomerville et al.l 
l2008al . l2012l : |Porter et al.ll20l4l . The GK and BR models 
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Figure 8. Same as Fig. [2] except here we compare the three “fiducial” models, GK, BR, and KS. We find remarkably similar 
results among all three models. The largest differences are in the predicted overall gas fraction and the H 2 fraction at high redshift 
in the lowest mass halos. 


are the same as the mo dels used in lPonning et alJ (|2014cJ l 
and lBerrv et al.1 d2014l l. Selected results for other model 
variants are shown in the Appendix. 

In Fig. |TT] we present predictions for the stellar mass 
function of galaxies from 2 = 0 to 2 = 6. We com¬ 
pare these predictions with a compilation of observa¬ 
tions as described in the figure caption. The similar¬ 
ity of the three model predictions is striking, particu¬ 
larly on the low mass end where we might have expected 
the differences to be largest. The most noticable differ¬ 
ences are instead at high masses, particularly at red- 
shifts z it 2. Here, the KS model produces significantly 
lower number densities of massive galaxies at 2 ^ 1, with 
the deviation growing with increasing redshift. It is im¬ 
portant to note that we have not accounted for the ex¬ 
pected errors in the observational estimates of the stel¬ 


lar masses when comparing to the models — this would 
tend to lead to an apparent increase in the numbe r of 
mass ive galaxies due to Eddington bias (see iLu et all 
120141 . e.g.). However, even if this effect were included, 
the earlier formation of massive galaxies predicted by the 
GK and BR models is clearly in better accord with re¬ 
cent observations. All three models suffer from the famil¬ 
iar excess of low-mass galaxies at 2 ~ 0.5-2 which, as 
we have already discussed, is a widespread problem in 
both semi-analytic models and numerical hydrodynamic 


simulations ( 

Fontanot et alJ 2009|; 1 Weinmann et al] 2012; 

White et al. 

2014: Somerville & Dave 20141. One of the 


important conclusions of this paper is that varying the 
star formation efficiency according to physically moti¬ 
vated recipes does not appear to he able to cure this 
problem within the current model framework. A possi- 
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Figure 9. Same as Fig. [2] except here we compare the results of our fiducial GK model with combinations of model ingredients 
that are similar to those used in several models in the literature. Models that neglect the effect of a varying UV background predict 
later star formation, higher cold gas fractions, lower H 2 fractions, and later chemical enrichment in low-mass halos. Metal-enhanced 
winds, when coupled with metallicity-dependent gas partition recipes, further delay star formation and enrichment. 


bly related problem is that over this redshift range, 
the predicted gas fractions in low mass galaxies in 
these same models may be too low dWh ite et al .ll2014l ; 
ISomerville fe Davi I20l4 iPopping et al.l l2014d . Popping 
et al. in prepj3- 

Fig. [12] shows a related quantity, the stellar frac¬ 
tion (stellar mass divided by halo mass; /star = m*/Mh) 
over the same redshift range. Our model predictions are 
now compared with constraints from ^sub)-halo abun¬ 
dance matching from lBehroozi et al.l (|2013i l. The conclu- 

1 An important caveat, however, is that the gas mass esti¬ 
mates that lead to this conclusion are based on indirect meth¬ 
ods. It remains to be seen whether this is conformed by direct 
observations of cold gas in these low-mass galaxies. 


sions are similar to the ones above, unsurprisingly since 
the /star constraints are derived from observational esti¬ 
mates of stellar mass functions (though not exactly the 
same ones plotted in our Fig. uni. The median stellar 
fractions are nearly identical in the three models in low 
mass halos (logM^/Mg ;$ 11) and are very similar in the 
GK and BR model over the whole range of halo masses. 
The median value of /star in the KS model is much lower 
in massive halos than in the other two models, and the 
difference increases with redshift to about 0.4-0.5 dex at 
2 = 6. 

Although the median values of /star are similar in 
the three models, the distributions differ significantly for 
low mass host halos. Fig. [13] and [14] show the distribu¬ 
tion of /star in halo mass and redshift bins, for central 
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Figure 10. Relation between total cold neutral gas density (Hi + H 2 ) and star formation rate density (SFRD). The solid black 
line shows the input H 2 -based star formation recipe (Big2). Slanted dotted lines show a star formation efficiency of 1%, 10%, and 
100% per 10 8 yr. Gray contou rs show observational estimates from the 13 Spiral galaxies from the THINGS+Heracles sample 
presented in lLerov et al.l ( 20081 '). Colored dots show a selection of 25 galaxies in our fiducial GK (left) and BR (right) models (at 
2 = 0), with a stellar mass range chosen to match the THINGS sample. Each point shows the value in an annulus with radius 500 
pc. The points are color-coded with the average gas phase metallicity of the galaxy. Note that in both models, galaxies with lower 
metallicity gas have lower SFRD for a given total gas density, because a smaller fraction of the gas is predicted to be in the form 
of H 2 . In the GK model, a direct dependence of H 2 fraction on metallicity is assumed. In the BR model, a dependence on disk 
midplane pressure is assumed, but this quantity turns out to be highly correlated with metallicity in our models. 


and satellite galaxies respectively. For all models and at 
all epochs, the distribution of /star becomes broader and 
more skewed with decreasing halo mass. For massive ha¬ 
los, the width of the distribution becomes slightly nar¬ 
rower with increasing time, while for low mass halos, a 
more noticable tail towards lower values of /star devel¬ 
ops with time. In the two intermediate halo mass bins, 
this tail is more prominent in the GK models than in 
the other two models. The predicted broadening in /star 
has potentially important implications for empirical halo- 
based models, which generally assume a narrow and fixed 
scatter in /star (Ah,). We show the results for these two 
types of galaxies separately because a) central and satel¬ 
lite galaxies are treated differently in SAMs. For example, 
in our models, satellite galaxies are not allowed to accrete 
new gas from the IGM; b) our predictions for differences 
between stellar fractions for satellites and centrals can 
provide useful input to empirical models such as Halo 
Occupation Distribution (HOD) and abundance match¬ 
ing models. Many such models do not distinguish between 
satellite and central galaxies. 


3.3.2 Star Formation rates and gas depletion times 

Fig. [T5] shows the specific star formation rate (sSFR 
= m,/m H tar) as a function of stellar mass over the red- 
shift range 2 = 0-6. Our model predictions are compared 
with a compilation of observations as described in the fig¬ 


ure caption. We have selected only “star forming” galax¬ 
ies using the criterion sSFR > 1/(3t.tr ( 2 )), where ttr(«) is 
the Hubble time at the galaxy’s redshift. This has been 
shown to produce similar results to commonly used obser¬ 
vational methods for selecting star forming galaxies (e.g. 
lLang et al.l 120141 ). Our models agree well with the ob¬ 
served slope and normalization of the Star Forming Main 
Sequence (SFMS) at 2 ~ 6-4, but then the normalization 
of the model SFMS drops below the observationally es¬ 
timated one between 2 ~ 3-0.5. At 2 ~ 0, the predicted 
SFMS has approximately the correct normalization for 
massive galaxies (here the precise value may be impacted 
by the details of the selection of “star forming” versus qui¬ 
escent galaxies), but the slope is much shallower than the 
observations suggest. This is anoth er facet of the “dwarf 
galaxy conundrum” discussed in I White et al.l d2014l ). and 
again is co mmon to most cosmologic al models of galaxy 
formation (ISomerville fe Pavel 120141 ). Our results show 
that this relation is extremely robust to changing the 
star formation recipe in models. 

Fig. [T7] shows the total gas depletion time td ep = 
(mm + m h 2 )/hr* in the fiducial GK and BR models, and 
in the GK+Bigl model. Here we compute the depletion 
time using only the SFR due to the ‘disc’ mode of SF, 
i.e. not including star formation due to merger-triggered 
bursts, but the plot looks very similar when the burst 
mode is included. Observations of nearby galaxies show 
that tdep increases with decreasing stellar mass, i.e. the 
conversion of cold gas into stars is less efficient in low 
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Figure 11. Stellar mass function evolution with redshift. Symbols show observational estimates as follows. In the z = 0.1, 1, 
2, and 3 panels, black square s ymbols show a doub le-Schecht er fit to a com pilation of observational e stimates. Observation s 
included in the fit are: z = 0.1 - iBaldrv et al.l (12008 h [Moustakas et all J2013T) ; z = 1 and z = 2 panels - iTomczak et all (|2014f) , 
iMuzzin et al.1 (|2013l ). z = 3 panels - iMuzzin et alJ ( 20131 ). The fits shown at z = 1, z = 2 and z = 3 are interpol ated to these 
redshifts from adjace nt redshift bins in th e original published results. In the z = 3 panel we als o show estimates fromlSantini et al.l 
(|2012|, triang l es) an d lCaputi et al.l d201 ll , crosses). In the z = 4 panel we sh ow estimates from [Duncan et al. ] <20l4 triangles) and 
I Caput i et alJ <|201ll . crosses). In the z = 6 panel we show the estimates from lDuncan et~al ] ll20li triangles). The purple solid line 
shows the results of the fiducial GK model, the orange dashed line shows the fiducial BR model, and the green dotted line shows 
the KS model. 


mass galaxies (e.g. iLerov et al.l 120081 ). A similar trend is 
indicated by the empirical estima tes of total gas depletion 
time from lPopping et al.l d2014al ). also shown in Fig. [17| 
for comparison. The empirical estimates are based on a 
SFR-halo mass relation inferred from abundance match¬ 
ing, and an indirect estimate of the Hi and H 2 masses 
from inverting the SFR density. These estimates rely on 
a number of assumptions (e.g., that disk cold gas radial 
profiles are well-represented by exponentials), and on the 
observed relationship betwe en size and stellar ma ss in 
disk-dominated galaxies. In IPopping et ail (l2014al ). the 
empirical predictions are shown only up to z ~ 3, because 
it is not known whether these assumptions and empiri¬ 
cal relations hold at higher redshift. Here we show the 
results of extrapolating the same method to 2 ~ 6, but 
these should be considered highly uncertain. 

Our three fiducial models reproduce the same qual¬ 
itative trends indicated by the observations and by the 
empirical predictions. First, depletion times are longer in 
lower-mass galaxies. In detail, the physics that is respon¬ 
sible for this trend is different in the three fiducial mod¬ 
els. Low mass galaxies have lower gas and stellar surface 
density on average. In the KS model, gas below the crit¬ 
ical surface density is not allowed to make stars, and low 
mass galaxies tend to have a larger fraction of their gas 


below this threshold. Low mass galaxies also tend to have 
lower gas-phase metallicities, and in the GK model, this 
results in less efficient formation of H 2 and thus of stars. 
In the BR model, the lower SFE in low-mass galaxies is 
due to their lower stellar surface density, which leads to a 
lower disc mid-plane pressure, and again a lower H 2 frac¬ 
tion. Second, in all models, tde P at a given stellar mass 
was lower in the past and increases with cosmic time. 
The GK and BR models show more pronounced evolu¬ 
tion and shorter depletion times (higher SFE) at high 
redshift, particularly in massive galaxies. This is due to 
the steeper dependence of the SFR density on gas den¬ 
sity adopted in our GK and BR models (slope IVsf = 2 
instead of 1.4). High-redshift galaxies contain higher sur¬ 
face density gas overall, and so the results are more sen¬ 
sitive to the slope of the SF relation at high gas densities. 
This result explains the less efficient formation of stars 
in massive galaxies at high redshift in the KS model rel¬ 
ative to the GK and BR models, seen in Fig. Illl and fTTl 
Note that already by z ~ 3, and increasingly so at higher 
redshift, the predictions using the Bigl recipe for SF are 
inconsistent with the empirical constraints. This suggests 
that the assumption of a constant H 2 depletion time in 
galactic disks (which is inherent in the Bigl recipe) may 
not be universally applicable. 
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Figure 12. The stellar mass of central galaxies divided by the total mass of their dark m atter halo, in redshift bins from z = 0-6. 
Dark gray lines and shaded areas show constraints from halo abundance matching from iBehroozi et ah f201 'l l. The purple solid 
line shows the results of the GK model, the orange dashed line shows the BR model, and the green dotted line shows the KS 
model. 


Fig. [18] shows the H 2 depletion time (tde P ,H 2 = 
mH 2 /m») in the GK, BR, and GK+Bigl models (the 
KS model is not shown, as we do not track H 2 self- 
consistently in this model). This figure, in combination 
with results shown in PST14, shows that the trends seen 
in our models in Fig. [TT] are due to a combination of two 
factors: a) at a given redshift, more massive galaxies have 
larger fractions of their cold gas in the form of H 2 , and at 
fixed mass, higher redshift galaxies have higher H 2 frac¬ 
tions; b) the H 2 depletion time is also shorter in more 
massive galaxies and at high redshift. 

3.3.3 Mass-metallicity relations 

Fig. [TU] and Fig. [20] show the mass-metallicity relation 
(MZR) for stellar and cold gas phase metallicities, re¬ 
spectively. Recall that the chemical yield parameter in 
our models has been adjusted to approximately repro¬ 
duce the n ormalizati on of the stellar MZR measured by 
iGallazzi et alJ (l2005lf . Our models naturally predict a 
slope for the stellar MZR that i s in fa i rly good agreement 
with observations (IWoo et al.l 120081 : iKirbv et al.l l2013f) 
down to very low stellar masses (m s tar ~ 10 7 Mq). The 
three fiducial models make very similar predictions for 
the stellar phase MZR, except that massive galaxies be¬ 
come enriched much earlier in our two new models (GK 
and BR) than in the KS model. This is owing to the 
steeper slope of our Big2 star formation recipe at high 
gas densities (see discussion in 1 )3.11 especially Fig. [6] and 


above). It is also interesting to note that our models pre¬ 
dict a much smaller dispersion in the stellar MZ R than 
the ob servational dispersion estimated bv IGallazzi et all 
(I2005T) . 

In Fig. 1201 the model results shown are for central, 
star forming galaxies selected using the same criteria de¬ 
scribed above. This is because the observational estimates 
of gas-phase metallicity are based on emission line diag¬ 
nostics from Hu regions, which are only detectable in 
star forming galaxies. A compilation of observational es¬ 
timates for the gas phase MZR is shown. We have con¬ 
verted the observed values of 12 + log(Q/H) to ZjZp, as- 
suming 12 + log(0/H) = 8.76 for the Sun llCaffau et al.l 
l201lll . and Z g jZQ = (O/H)/(O/H) 0 . Note that the ob¬ 
served MZR is uncertain by a factor of 2-3 as different 
calibration methods produce d ifferent zero-points and , 
to some extent, different slopes (IKewlev fe Ellisonll2008i l. 
The predictions of our three models are again quite sim¬ 
ilar, except that the KS model again produces later en¬ 
richment of massive galaxies, so the MZR is shallower 
at 2 ~ 3-6. All three models produce a cold gas phase 
MZR that is, taken at face value, considerably steeper 
than the observed gas phase MZR. Moreover, the pre¬ 
dicted evolution of gas phase metallicity in our models 
is quite different from that implied by current observa¬ 
tions. The models predict that the gas phase metallicity 
for galaxies of fixed stellar mass declines slightly with de¬ 
creasing redshift, while observations indicate an increase 
of almost a factor of two between 2 ~ 2.2 and 2 ~ 0. This 
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Figure 13. Distribution functions for the stellar fraction (/star = mstar/Affc) of central galaxies in bins of halo mass (log A-7), = 
9.25 — 9.75, 9.75 — 10.25, 10.25 — 10.75, 10.75 — 11.25) and redshift as indicated on the panels. The purple solid line shows the 
results of the GK model, the orange dashed line shows the BR model, and the green dotted line shows the KS model. 


discrepancy was shown previously by I White et alJ J2014I ) 
for our KS models; we see here that the results are qual¬ 
itatively similar for our new fiducial GK and BR models. 
We discuss possible reasons that our models reproduce 
the stellar MZR fairly well but seem to fail to reproduce 
the observed gas phase MZR in H4.2I 


3.3.4 Evolution of Global Quantities 

In Fig. [21] we show our model predictions for the evolu¬ 
tion of the global SFR density, global stellar mass density, 
and average metallicity of cold gas and stars over cosmic 
time. The figure shows that the three models predict al¬ 
most identical global SFR densities at low redshift, while 
at z it 2, the BR model produces the highest SFR density, 
and the KS model the lowest, with the GK model in be¬ 
tween. We also see that our model predictions are in rea¬ 
sonably good agreement with observations at z S 2 and 


6, though this is in part a fortuitous cancellation 
the models overproduce galaxies with low SFR but under¬ 
produce ones with high SFR. Moreo ver, the observational 
estima tes of the SFR density from iMadau fe Dickinson! 
(120141 ) are integrated only down to 0.03 L„, while our 
theoretical predictions are integrated over all galaxies. 
This same behavior is echoed in the build-up of the global 
stellar mass density. The largest difference in the three 
models appears in the evolution of the stellar and cold 
gas phase metallicity. The models differ in the normal¬ 
ization and evolution of the mean metallicity for gas and 
stars. In addition, in the KS model, the mean metallicity 
of cold gas and stars is very similar, while there is a much 
larger difference between the stellar and gas phase metal- 
licities of gas and stars in the BR model and GK models. 
As one can see in Fig. 1191 and 1201 at stellar masses above 
7Tistar ~ IQ 7 Mq, the models make similar predictions for 
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Figure 14. The same as Fig. 1131 but for satellite galaxies (sub-halos). The halo mass here is the mass of the halo when it first 
becomes a sub-halo. 


the relative metallicities of gas and stars. The differences 
seen in Fig. [5T]are entirely due to very low mass halos. 


4 DISCUSSION 

4.1 Interpreting our results: the equilibrium 
model 

One of the main conclusions of our work is that modi¬ 
fying the recipes for how cold gas is converted into stars 
has very little effect on the properties of low-mass galax¬ 
ies (m* ^ M c har, where M c h ar is the “knee” in the stel¬ 
lar mass function). Instead, modifying the star formation 
recipe mainly changes the ratio of cold gas to stars in 
galax ies. Sim ilar con clusions were reached in the study 
by [White et alj (120141 ), in which more extreme (though 
in some cases less physically motivated) modifications to 
the SF recipe in similar models were made. This is be¬ 
cause in our models, star formation in low-mass galaxies 


is strongly self-regulated: if star formation is made less 
efficient, less gas is ejected by stellar winds, leading to 
more efficient star formation, and vice versa. 

A number of recent works have pointed out this 
property of self-regulation, which is a rather generic 
feature of modern galaxy formation models, arising 
from the broadly adopted hypothesis of efficient stellar- 
driven feedback ( Scha 2 e_et_al ] 201h| Haas et ah! [2013 :: 
IWhite et al .1120141 : Somerville fe Paval2014l ). A useful an¬ 
alytic framework for understanding the behavior of the 
fairly complex intertwined suite of physical processes at 
play in state-of-the-art galaxy formation simulations is 
the “equilibrium model”, sometime s called the “bath¬ 
tub model” (e.g. I Paveet al.1 120121 : iDekel et al.1 l201.il : 
iDekel fe Mandelkerl 12014 ). The basic assumption in this 
model is that due to self-regulation, on some timescale 
f eq , galaxies establish an equilibrium state in which the 
rate of change of their cold gas reservoir is small, i.e. 
m C oid — 0. Once equilibrium is established, the star for¬ 
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Figure 15. Mean specific star formation rate (sSFR = ra*/m*), as a function of stellar mass for our three fiducial models (purple 
solid: GK; orange dashed: BR; green dotted: KS), in redshift bins from z = 0-6. The blue contours show the conditional sSFR in 
the GK model. The horizontal gray line shows the sSFR corresponding to 1/(3 £#), where tn is the Hubble time at that redshift. 
Only galaxies with s SFR > 1 /(3 tn) are included in the mean. Symbols show a compilation of observations for star forming galaxies 
as follows: z = 0.1 - Salim et al. ( 20071. open circ les); z = 1 - 1 Whitaker et ah (|2014| . pentagons, i nterpolated in redshift from the 
published results) ; z = 4- [Steinhardt ~et al.l <2014 crosses); * = 4 and z = 6 - 1 Salmon et alJ <20l4 circles); all panels - fit to data 
compilation from [Soeagle et~ J <2014 squares). 


mation rate is balanced by global inflows and outflows, 
rh t = rhin/(l + ??), where rhi n is the rate at which gas 
flows into the galaxy due to cosmological accretion and 
rj = rhout/rh, is the mass loading factor of a large-scale 
stellar driven outflow. The time for a galaxy to come into 
equilibrium (or to re-establish equilibrium after a disrup¬ 
tion) is 


where tsf = m C oid /m* is t he star f ormation (or gas de¬ 
pletion) timescale lIDave et al .112013 . hereafter DF012). 

In simulations, r/ is a fairly strong inverse function 
of galaxy mass, while tsf is a weaker function of galaxy 
mass (DF012). Therefore low mass galaxies come into 
equilibrium earlier. This helps us to understand why 
changing our star formation recipe had less impact on 
low-mass vs. massive galaxies. In addition, it explains 
why high-mass galaxies were affected only at high red- 
shifts - at z it 2, these galaxies have not yet come into 
equilibrium. This work makes the interesting prediction 
that observations of massive galaxies at very high red¬ 
shift (z it 2) will place strong constraints on the physics 
of star formation in cold dense gas, while constraints on 


the low-mass end of the galaxy stellar mass function at 
high-z will mainly constrain the physics of outflow^. 

4.2 Mass-metallicity relations for gas and stars 

It is puzzling that our models reproduce the stellar MZR 
fairly well but predict a much steeper gas phase MZR 
than observations appear to indicate. This has been seen 
in other models as well — models that invoke a weaker de¬ 
pendence of mass outflow rate on galaxy circular velocity, 
and normalize their yield parameter to the observed gas 
phase MZR, produce better agreement with the observed 
gas p hase MZR but then fail to reproduce the stellar 
MZR dLu et al.ll2014h . |Peeples fe Somervillel J2013I ) com¬ 
bined empirical star formation histories derived from the 
observed SFMS with the observed relation between SFR, 

2 All of this discussion implicitly assumes that “preventative” 
feedback — physical processes that could prevent gas from ac¬ 
creting into galaxies or becoming available for star formation 
— is sub-dominant. At very low masses, preventative feedback 
due to photo-ionization squelching likely becomes important. 
Preventative feedback due to AGN heating and winds, and 
virial shock heating, is probably dominan t in massive galax¬ 
ies at late times (z t 2). See DF012 and ISomerville &; Davd 
<20l4) for a more complete discussion. 
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Figure 16. The conditional probability distributions of sSFR in stellar mass bins at different redshifts, for our three fiducial 
models (purple solid: GK; orange dashed: BR; green dotted: KS). 


stellar mass and gas phase metallicity (assumed to be uni¬ 
versal) and predicted the stellar MZR, finding fairly good 
agreem ent with observations. However. iMunoz fe Peeplea 
d2014h did a more realistic calculation using a similar ap¬ 
proach, but accounting for stochasticity in the SF histo¬ 
ries and quenching, and found more significant tension 
between the gas and stellar phase MZR. 

We can see by comparing Fig. [T9] and [20] that our 
models predict that the stellar metallicity in galaxies is 
about a factor of 1.5-1.7 lower than the gas phase metal¬ 
licity, with weak trends on stellar mass and redshift. We 
can also see from the figures in m\ that the gas phase 
metallicity tends to increase rapidly and monotonically 
with time in our models. The stellar metallicity is ef¬ 
fectively a mass-weighted average over the chemical en¬ 
richment history of the galaxy, so it makes sense that the 
stellar metallicities are slightly, but not enormously, lower 
than the gas phase metallicities at any given time. Taken 
at face value, the observational results — which imply 


that Zg as /Z s t a r is as high as a factor of ~ 10 or more, 
and is a fairly strong function of stellar mass — may 
be difficult to reproduce in cosmological models without 
invoking accretion of highly metal pre-enriched gas. 

An alternative explanation is that the normalization, 
and possibly the slope, of the gas and stellar phase MZRs 
are not accurately calibrated to the same system. Indeed, 
some gas phase metallicity indicators do yield a gas MZR 
normalization and slope that is more c onsistent with the 
stellar MZR (jKewlev fc Ellisonl 120081 : IMunoz fc Peeplesl 
120141 ). Another potential issue is that we have plotted 
stellar mass weighted metallicities, while the observed 
stellar metallicities are luminosity weighted. However, 
other studies have found that the luminosity weighted 
stellar MZR does not differ signi ficantl y in slope fro m the 
stell ar mass weighted MZR ( Trager fc Somervillell200 9i: 
IPeeples fc Somerville!l2013! 'l. 

Another issue is that different observational probes 
are sensitive to different chemicaf elements. The stellar 
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Figure 17. Total gas depletion time, defined as the total cold neutral gas mass (Hi + H 2 ) divided by the star formation rate. 
The purple solid lines show the predictions of the fiducial GK model, the orange dashed lines show the BR model, and the green 
dotted lines show the KS model. The lavender dot-dashed lines show the GK+Bigl model. The 16th, 50th, and 84th percentiles 
are shown for central star forming galaxies in the models. The horizontal black dashed line shows the age of the Universe at that 
redshift . Th e open diamonds in the z = 0 panel show observational estimates for galaxies in the THINGS +Heracles sample f rom 
I Leroy et al.l (120081) . The horizontal gray line shows the average molecular gas depletion time estimated bv lLerov et alj <120131) for 
nearby galaxies; the shaded gray area indicates the uncertainty in the measurement due to the uncertain conver sion factor between 
CO and H 2 . The gray circles show the estimates obtained via the empirical method of lPopping et all d2014al) . These are plotted 
with open symbols at z > 3 to indicate that the estimates are quite speculative in this redshift regime (see text for more details). 
All three models reproduce the observed trend of decreasing depletion time with increasing stellar mass, but different underlying 
physics are responsible for the trends in the different models. 


metallicities derived by iGallazzi et al.l (|20051 ) are sensi¬ 
tive to a combin ation of Fe a nd M g, and the stellar 
MZR derived by Kirby et al.l (|2013l ) measures [Fe/H]. 
IGallazzi et al.l (120051 1 quote their results in terms of 
•Zstar/Z© and claim that their results are independent of 
a/Fe (A. Gallazzi, priv. comm.). The observational gas 
phase MZRs are primarily sensitive to a elements and are 
quoted in terms of oxygen abundance (12 + log(0/H)). 
We have assumed that Z gas /Zo = ( Q/H)/(Q/H ) p) and 
Zst&rlZ q = (Fe/H)/(Fe/H) 0 for the IlGrby et all J2013fl 
observations. This is equivalent to assuming that all the 
stars in our model galaxies have (a/Fe) = (a/Fe)©. How¬ 
ever, (a/Fe) is known to differ significantly fromJheSolar 
value in stars in our own Galaxy fe.g. ISto ll et al] l2013h . 
in nearby dwarf galaxies JToIstov et al.l 120091. and ref¬ 
erences therein), and in giant ellipticals (|Thomas et all 
120051: llrager et alJl200d) . 


In the models presented here, we track the total 
metallicity assuming a constant yield, and we also adopt 
the instantaneous recycling approximation. Metals are 
produced in direct proportion to the formation of stars, 


so enrichment in our models probably most closely traces 
a elements, but we normalized our yield parameter to 
observations that are also sensitive to Fe (see above). 
This simple version of single-element chemical enrich¬ 
ment with the instantaneous recycling approximation is 
the standard approach adopted in semi-analytic models. 
However, physical processes in the models are actually 
dependent on different elements in potentially significant 
ways. For exampl e, the (also widely ad o pted in SAMs) 
cooling tables of ISutherland fc Dopital dl993h used to 
model the cooling rate of hot halo gas are parameter¬ 
ized via [Fe/H] and adopt assumed [Fe/H]-dependent 
abundance ratios. However, H 2 formation is more closely 
tied to a elements such as oxygen and carbon, which 
are primary coolants in the interstellar medium (e.g. 
iGlover fe C lark 2014). One conclusion of the work pre¬ 
sented here is that if we wish to include more realis¬ 
tic physics in our models, it is important to track mul¬ 
tiple chemical elements and their production via dif¬ 
ferent channels (stars of various masses, Type II SNae 
and prompt and delayed Type la SNae) and on dif- 
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Figure 18. H 2 depletion time, defined as the H 2 mass divided by the star formation rate. Models shown (colored lines) are as 
in Fig. in The horizontal black dashed line shows the age of the Univ erse at that redshift. The horizontal gray line shows the 
average molecular gas depletion time estimated by iLerov et al.l ll2013h for nearby galaxies; the shaded gray area indicates the 
uncertainty in the measurement due to the uncertain conversion factor between CO and H 2 . Our models predict that at high 
redshift, molecular gas depletion times were significantly shorter, and there is a much stronger trend between galaxy stellar mass 
and H 2 depletion time. 


ferent timescales. Several groups have developed SAMs 
that include more detailed chemical e volution models 
llArrigoni et al.l l20ld : lYates et al.l l2013h . We have inte¬ 
grated the more soph i sticat ed chemical evolution models 
presented in lArrigoni et al] (l2010h within our new SAMs, 
including the new metallicity-dependent H 2 formation 
and H 2 -based star formation recipes as described here, 
and plan to investigate the predicted metallicities of cold 
gas and stars and their evolution in these models in a 
future work (Peeples et al. in prep). 

In the meantime, we can perform an empirical cor¬ 
rection for the variation of a/Fe to see if this is a plausi¬ 
ble explanation for the discrepancy. If we assume that 
the metallicity tracked in our models is a ctua lly Fe, 
and apply the empirical relation presented bv lstoll et al.l 
(12013r n to our model galaxies to “convert” to [O/H], we 
find much better agreement between our predicted gas 
phase MZR and at least some calibrations of the ob- 
ser ved MZR (see Fig.|20ll. N ote that, conceptually follow¬ 
ing iMunoz fe PeeplesM20141 ) , we are effectively assuming 
that a relationship between [Fe/H] and [O/H] derived for 
individual stars in the Milky Way holds for the average 
stellar and gas metallicities in galaxies with a variety of 


3 [Fe/H] = -0.34 + 1.25[0/H] 


star formation histories — an assumption that may well 
not be valid. However, it suggests that properly account¬ 
ing for non-Solar a /Fe and its possible trends with other 
galaxy properties (such as stellar mass and metallicity) 
may at least partially relax the tension between the stel¬ 
lar and gas phase MZR seen in our models and others. 


4.3 Caveats and limitations of our models 

Understanding how the “small scale” processes of star 
formation and stellar feedback interact with cosmological 
scale processes such as galactic scale inflows and outflows 
to shape the observable properties of galaxies is currently 
one of the major unsolved problems in the study of galaxy 
formation and evolution. The models presented here ne¬ 
glect a large number of processes that are thought to 
be important in influencing how efficiently gas can form 
molecules, and in turn how stars form within molecular 
gas. For example, we do not consider the possible impact 
of the local shear field, and non-axisymmetric perturba¬ 
tions such as spiral arms and bars. Nor do we attempt 
to model the “local” effects of stellar feedback (through 
stellar winds, supernovae, and Hn regions) on star forma¬ 
tion. We instead assume that the efficiency of converting 
molecular gas into stars is roughly constant, as suggested 
by observations of nearby spiral galaxies (IBigiel et al.l 
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Figure 19. St ellar mass vs. stellar metallicity. Five-pointed star symbols and dashed lines show observational estimates from 
Gallazzi et al.1 (l2QQ5h . and six-pointed stars show the fit to the observed MZR estimated from Local Group dwarf galaxies by 
Kirb y et a l. tom The purple lines show the 16, 50, and 84th percentiles for our fiducial GK model, the orange line shows the 
BR model, and the green line shows the KS model. 


l2008i. l201ll) . As pointed out by iKrumholz et al.l d2012l i 
and many others, the efficiency of forming stars within 
GMCs is surprisingly low, only about 1% per free fall 
time. Our picture is that local feedback processes are re¬ 
sponsible for setting that efficiency, and that when we 
smooth over several 100 pc regions of the ISM, it aver¬ 
ages out to a nearly universal value. 

However, there have been some recent studies that 
suggest that the galaxy-averaged value of this molecular 
star formation efficiency (often expressed as a depletion 
time, fde P ,H 2 = mjr 2 /rh») may vary significantly from 
galaxy to galaxy, and may ha ve a strong d epe ndence on 
global galaxy properties. iSaintonge et al .1 (1201111 find that 
in the COLD GASS sample, tde P ,H 2 is weakly correlated 
with the galaxy stellar mass and stellar surface density, 
and rath er strongly (anti-) correlated with sSFR. How¬ 
ever, ISaffitonge^Tahj (1201 il l adopted a constant (Galac¬ 
tic) value for the conversi on factor between CO and H 2 
(qco). ILerov et all (hoillh found similar correlations in 
their sample of 30 disk galaxies from the HERACLES sur¬ 
vey, but found that most of the correlation disappeared 
when they applied a theoretically motivated dependence 
of qco on the dust-to-gas ratio. They found smaller resid¬ 
ual variations in tde P ,H 2 , mainly associated with nuclear 
gas concentrations. This suggests that there may be dif¬ 
ferent values of fd ep ,H 2 in undisturbed disk galaxies and in 
galaxies experiencing mergers and in teractions (see also 
iDaddi et al.l l2Qlol : iGenzel et al.ll2010ll . possibly due to a 
super-linear dependence of SFRD on H 2 density, as as¬ 


sumed in our fiducial models. We do include an enhance¬ 
ment of SFE in mergers in our models, but the treatment 
is based on hydrodynamic simulations of binary mergers 
with a rather outdated treatment of sub-grid physics, so 
this is clearly an area that should be explored with state- 
of-the-art high-resolution hydrodynamic simulations. 

Another limitation of our approach is that although 
we compute the H 2 fraction and Esfr in radial annuli 
in each disk, and integrate over the disk to obtain the 
global properties, we do not store the information on 
the stellar, gas, and metal content of each annulus over 
different timesteps. Therefore we assume that both the 
gaseous and stellar disks have radial exponential profiles, 
and adopt a simple fixed factor relating the size of the 
gaseous and stellar disk. While this may be a reasonable 
approximation on average, it may miss important trends. 
We also use the global values of the gas phase metallicity 
and SFR (which we use as a proxy for the UV radiation 
field) in the GK recipe, instead of the local values of these 
quantities in annuli. In future work, we plan to construct 
more detailed models of disks in which all of these quan¬ 
tities are t racked as a func tion of radius, along the lines 
of work by iFu et all (l2010l l and iDutton et all (120101 1 . 

4.4 Comparison with previous work 

Several othe r group s hav e carr ied out studies similar 
to this one. lLagos et al.l (j2011bl l considered two models 
without gas partitioning. The first unpartitioned model 
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Figure 20. Stellar mass vs. cold gas phase metallicity. B l ack sy mbols show observational estimates of the gas-phase metallicity: 
z = 0.1 - 1 Peeples et al.ll (1 20141 . p luses)• [ Andre ws &; M artini (2 013 . filled circles). In all panels, the filled squares show the compilation 
of observational estimates from Zahid et al. 1 120131) . For comparison, we also show the observational estimates of stellar metallicity , 
as in Fig. m with gray star symbols. Colored lines show model predictions, as in Fig. 1191 The dot-dashed purple line shows the 
GK model with an approximate correction for varying [a/Fe] (see text). Taken at face value, our predicted gas-phase MZR appears 
to be much steeper than the observational estimates. However, properly accounting for varying abundance ratios of a versus Fe 
elements may at least partially remove this tension (see text). 





Figure 21. Global history of star formation, stellar mass assembly, and metallicity as a function of redshift. The solid purple 
line shows the results of our fiducial GK model, the dashed orange line shows the BR model, and the dotted green line shows 
the KS model. Left: global star formation rate density; gray lines show the fit to the compilation of observational estimates from 
iMadau &; Dicki nsonl (12 0141). Middle pa nel: global stellar mass density; symbols show selected observational estimates taken from 
Table 2 of lMadau &; Dickinsonl l)2014l ). Right: mean mass-weighted metallicity; thick lines show the metallicity of the cold gas 
component and thin lines show the stellar component. The diamond symbol at z ~ 0 is the mean stellar metallicity derived by 
lOallazzi et al.l d2005l) . 


adopts the origina l SF relations implem ented in the 
iBaueh et al.l d2005l ) and iBower et al.l (120061 ) GALFORM 
models, in which the SFR was assumed to be proportional 
to the tot al cold gas mass divided by a timescale r*. In the 
iBaugh et al.l l|2005l ) models, r* was scaled with the galaxy 


circul ar velocity to a power, and in the IBower et al.l 
(120061 ) models, also with the galaxy dynamical time. The 
second unpartitioned model used a Kennicutt-Schmidt 
SF recipe similar to the one we have adopted here. They 
also considered two models with gas partitioning: one 
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with the pressure-based BR recipe, and one with the 
metallicity-based KMT recipe. In their BR model, they 
adopted SF relations similar to our Bigl and Big2 recipes. 
In the KMT model, they used the SF relation given by 
KMT. They did not attempt to separate the effects of 
the different gas partitioning recipes versus H 2 -bas ed SF 
recipes. We focus on the results of the lBower et al.l |200fl) 
variant of the GALFORM models, which are more sim- 
ilar to our models t han the iBaugh et al.l d200fif l variant. 
lLagos et al.i ( 2011b! ') do not show stellar mass function 
predictions, but find that the rest frame K-band luminos¬ 
ity function at 2 = 0, 1, and 2 shows little change between 
the six different star form ation recipes they explored, con¬ 
sistent with our results. lLagos et al.l d2011bl l emphasize 
differences in the SFR distributions as a function of stel¬ 
lar mass, in particular the prominence of a passive popu¬ 
lation in models with different star formation recipes. We 
find very small differences in the SFR distributions be¬ 
tween our different models, and note that the prominence 
and location of a passive population will certainly also be 
very sensitive to the trea tment of AG N feedba ck. It is also 
interesting to note that lLagos et alJ (! 201 lal ') end up fa¬ 
voring their BR models and strongly disfavor the KMT 
recipe because they find that it does not reproduce the 
observed Hi and H 2 gas scaling relations. However, we 
showed in PST14 that both our BR and GK models do 
about equally well at reproducin g available observat ions 
of cold neutral gas in galaxies. lLagos et al.l (2011b! ') do 
not sh ow me tallicity predictions. 

iFu et all <201211 conducted a similar study, investi¬ 
gating four SF recipes and separately studying two gas 
partitioning recipes, BR an d KMT, within t he framework 
of the models developed bv Fu_et_alJ fl2010|) based on the 
MPA Millennium SAMs Icroton et al.1 200fil : lOuo et all 
l201ll ') . Their “Bigiel” SF recipe is similar to our Bigl 
recipe, but depends on the H 2 fraction, steepening below 
a critical value. Their “Kennicutt” recipe is similar to our 
KS recipe. Their “Genzel” recipe contains a linear scaling 
of Espr with E_h 2 , but also scales as the inverse galaxy 
dynamical time (which is redshift dependent). They also 
consider the KMT SF recipe. Again their results are quite 
consistent with ours. The stellar and H 2 mass functions 
are almost identical for all models, while the largest dif¬ 
ference between models is in the Hi content. This is the 
same conclusion that we reach based on PST14 and this 
study. All of their models underproduce massive galaxies 
at high redshift (z it 2), as we found with similar_S!F pre¬ 
scriptions to the ones they adopted. IFu et all (120121 ') do 
not show their m s tar-SFR relation, but they show that 
the different star formation recipes can lead to substan¬ 
tial deviation between model predictions for the cosmic 
SFR density at high redshift (z ji, 3-4), as we also find. 
Interestingly, IFu et all <2012i ) find that their gas phase 
MZR is too shallow compared with observations — the 
opposite problem to the one we encounter in our models. 
This is probably due to their adopted scaling for stellar 
driven winds. They assume that the mass outflow rate 
rhout oc m*, i.e., a fixed mass loading factor, while we 
assume that the mass loading factor rhout/m* scales ap¬ 
proximately with inverse circular velocity squared. This 
dependence of MZR slope on wind sca ling parameters 
is well known dPeeples fe Shankaj[20 11). They also find 


that the redshift evolution of the MZR is quite sensitive 
to the SF recipe adopted; in particular, recipes in which 
SFE scales with galaxy dynamical time predict very weak 
or no evolution in the MZR, while their models without 
dynamical time scalings predict stronger evolution. None 
of our models contain an explicit scaling with dynamical 
time, but the non-linear slope of our Big2 and KS recipes 
has a similar effect, consistent with the weak evolution in 
the MZR seen in our m o dels. 

iKrumholz fe Dekell d2012l ) implemented an updated 
version of the KMT recipe within simplified semi-analytic 
models that only follow the mass accretion history of the 
main branch, and do not track the full merger trees. In 
their fiducial model they additionally assume a constant 
mass loading factor for stellar-driven winds (no depen¬ 
dence on galaxy mass or circular velocity) and strongly 
metal-enhanced winds. When we implement similar in¬ 
gredients in our models (except that we retain our “en¬ 
ergy driven” stellar wind scalings), we obtain qualita¬ 
tively similar results. Namely, a metallicity-dependent 
formation efficiency for H 2 and self-consistent H 2 -based 
star formation recipe can significantly suppress and de¬ 
lay star formation and metal enrichment in very low-mass 
halos. We find that this effect is considerably stronger in 
our KMT+MEW models than in our fiducial GK models. 
This is because a) in our fiducial GK model, the effect of 
a varying UV background partially compensates for the 
metallicity dependence of H 2 -formation; and b) the metal 
enhanced winds delay enrichment of the cold gas, further 
suppressing star formation (see However, we stress 

that a noticable effect is seen only in halos with virial 
mass Mh S IQ 1 0 5 Mq and at early times (z it 1). There¬ 
fore, although IKrumholz fe Dekell (120121 ') do not show 
their predicted stellar mass functions or stellar fractions, 
it is likely that their model also still suffers from the 
overprediction of low-mass galaxies (m s tar ~ 10 9 ~ 10 M@) 
at intermediate redshifts (0.5 S z S 4) that we have dis¬ 
cussed extensively above (see also Fig. 1221 and 1231) . Cer¬ 
tainly they see the same qualitative problems with pre¬ 
dicted specific star formation rates being too low at in¬ 
termediate redshifts that we have described here. 

Galaxies hosted by halos in the strongly affected 
mass range have typical stellar masses of m s t ar ~ 
10 ,_8 M@, and it is probably not feasible to obtain com¬ 
plete samples of these objects at high redshift with ex¬ 
isting facilities. However, this strong suppression of star 
formation in low mass halos could have implications for 
reionization, future observations with the James Webb 
Space Telescope, a nd stellar archaeolo gy in local dwarf 
galaxies. Moreover, iBerrv et al.l d2014l ') showed that the 
strong suppression of H 2 formation and star formation in 
low-mass halos predicted by the KMT-like models would 
lead to a large population of “barren” halos that never ex¬ 
perience si gnificant star forma t ion, and so are fille d with 
Hi (see also Kuhlen et alj20l3 ). IBerrv et al.l (120141 ') found 
that the existence of this population is in apparent ten¬ 
sion with observations of Hi absorption systems at high 
redshift. 

We found that the inclusion of metal-enhanced winds 
(MEW) produces a steeper MZR, leading to an even 
g reate r tension with observations. Krumh olz fc Dekel 
(12012! ) also found that their models produce a steeper 
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gas phase MZR than observations, in spite of their 
adopted constant mass loading factor which gener- 
ally leads to shallower predicted MZR. Interestingly, 
iKrumholz fc Dekell ll20121 ) find significant evolution in the 
gas phase MZR from z ~ 2-0 in their models, in better 
apparent agreement with observations. 


5 SUMMARY AND CONCLUSIONS 

We have presented new models of galaxy formation, set 
in the framework of cosmological merger trees, which 
include physically motivated recipes for the partion- 
ing of cold gas into an atomic, molecular, and ionized 
phase. These models have several advantages: first, we 
can make explicit predictions for the atomic and molec¬ 
ular gas properties of galaxies over cosmic time, which 
can be directly confronted with observations from cur¬ 
rent and upcoming facilities. The first of these predic¬ 
tions, for Hi an d H 2 gas o bserve d in emission, were pre¬ 
sented in IPopping et al.l ll2014d) , and predic tions for Hi 
gas ob s erved in absorp t ion we re presented in lBerrv et al.l 
(12014h . [Popping et al.l ll2014bl ) extended these models to 
predict sub-mm line emission from several atomic and 
molecular species. These predictions may be directly 
compared with observations from current and upcoming 
facilities such as ALMA (Popping et al. in prep), and 
also used to plan future observations with these facilities. 
Second, we can implement more physically motivated H 2 - 
based recipes for star formation within our models. In this 
paper we focussed on the predictions for the stellar mass 
content, star formation rate, and metallicity of galaxies 
in our new models, and compared these with available 
observational estimates from z ~ 0-6. 

We summarize our main conclusions below: 

• We performed a number of tests of the robustness 
of our models to resolution and various parameter values 
and ingredients. We find that our code is very well con¬ 
verged with respect to variation of the mass resolution of 
our merger trees. In addition, our results do not depend 
sensitively on the assumed values of the metallicity of 
pre-enriched gas or the molecular hydrogen floor (within 
a reasonable range of values). 

• Accounting for a reservoir of ionized gas (due to an 
internal and external photo-ionizing radiation field) in 
our models, which is not allowed to form molecular hy¬ 
drogen or stars, does not significantly change our predic¬ 
tions. 

• We explored the effect of adopting different H 2 -based 
star formation recipes. All the recipes we considered 
gave similar results for low-mass galaxies. Models that 
adopted a “two-slope” recipe (Big2), which has a linear 
dependence of star formation rate density on molecular 
gas surface density below a critical value, steepening to 
Esfr oc Eh2 2 above a critical value of Eh2, produced 
more efficient star formation and metal enrichment in 
massive galaxies at high redshift. Models that implement 
the Big2 recipe appear to be in better agreement with 
current observational estimates of the number density of 
massive galaxies at high redshift. 

• We explored the effect of different recipes for par- 


tioning gas into an atomic and molecular component. The 
metallicity-based GK recipe and the pressure-based BR 
recipe gave surprisingly similar results, perhaps because 
of the strong correlation between disk mid-plane pres¬ 
sure (oc E*, to first order) and metallicity in our models. 
The KMT and GKFUV recipes, which do not include 
a dependence on the FUV radiation background as our 
fiducial GK recipe does, predicted less efficient forma¬ 
tion of H 2 , less star formation and metal enrichment at 
early times, and later stellar mass assembly. These differ¬ 
ences are only noticable, however, in very low mass halos 
(log(Af h /M 0 ) < 10.5). 

• Both of our new fiducial models (GK and BR) re¬ 
produce the curvature in the relationship between total 
cold gas surface density and Esfr seen in observations of 
nearby spiral galaxies. Our results illustrate the difficulty 
of disentangling the physical processes that are responsi¬ 
ble for the scatter in this relationship, however, because 
of the strong correlations in galaxy properties. 

• The stellar mass function and mean stellar fractions 
(/star = m*/Mh) of galaxies in our three fiducial models 
(the “classic” KS model, the metallicity-based GK model 
and the pressure-based BR model) are almost identical 
for low-mass galaxies at all redshift z ~ 0-6. Both of the 
new models (GK and BR) predict earlier formation of 
massive galaxies, in better agreement with observational 
estimates than models that adopt the KS recipe. 

• Although the median values of /star are very similar 
in low-mass halos in all three models, the models can 
have significantly differently shaped distribution func¬ 
tions P (/star |-M/). In particular, the GK model tends 
to have a much more pronounced tail to low values of 
/star in low-mass halos (log(M/,/MQ) £ 10.5). The pre¬ 
dicted broadening in /star has potentially important im¬ 
plications for halo occupation models, which generally 
assume a narrow and fixed scatter in /star(M^). 

• All three fiducial models produce nearly identical 
predictions for the relationship between stellar mass and 
SFR at all redshifts. Even the distributions of sSFR at 
a given stellar mass are very similar. The KS model pre¬ 
dicts a slightly narrower distribution of sSFR in high- 
mass galaxies at z X) 1 than the other two models. 

• All three models predict a weak dependence of the 
gas depletion time (tdep = (mm +mH 2 )/m.) on stellar 
mass, in agreement with observations of nearby normal 
disk galaxies and empirical estimates. The predicted fde P 
decreases by about 1.3 dex from 2 ~ 0 to z ~ 6. The KS 
model predicts milder evolution in the depletion time for 
massive galaxies to high redshift, resulting in longer de¬ 
pletion times in massive high redshift galaxies compared 
to the other two models. 

• All three models predict quite good agreement with 
the observed z = 0 stellar mass versus metallicity relation 
(MZR) for stellar metallicities, but predict a gas phase 
MZR that is much steeper than observational estimates 
taken at face value. However, this tension may perhaps 
be relieved by properly accounting for the possible de¬ 
pendence of [a/Fe] on galaxy properties. The KS model 
predicts later metal enrichment of massive galaxies, lead¬ 
ing to a shallower MZR at high redshift. In tension with 
observational results, all of our models predict a nearly 
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constant or slightly declining metallicity for galaxies se¬ 
lected at fixed stellar mass from z ~ 4-0. 
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APPENDIX: RESULTS FOR ADDITIONAL 
MODEL VARIANTS 

In this Appendix we show results for the stellar mass 
functions and stellar fractions in several additional model 
variants, in order to aid the interpretation of the results 
presented in the main text. In Fig. 1221 we show the stel¬ 
lar mass function at a = 0, 1, 2, and 6 for our fiducial 
GK model (the same one shown in Fig. 1111) . compared 
with the GK model with a fixed value of the UV radi¬ 
ation background Umw = 1 (GKFUV), the GK model 
for gas partioning with the Bigl SF relation (GK+Bigl), 
and a model with the KMT recipes for gas partioning 
and a KMT SF relation (see Table [5J. Fig. 1231 shows the 
median stellar fraction as a function of halo mass for cen¬ 
tral galaxies, in the same suite of models, with the addi¬ 
tion of the GK model with photo-ionization “squelching” 
switched off shown in the 2 = 0 panel only. 

These plots illustrate several points, which we al¬ 
ready discussed in 43.11 First, the metallicity and UV ra¬ 
diation field dependence in the GK recipe partially coun¬ 
teract each other (lower mass galaxies have lower metal¬ 
licity, resulting in less efficient H 2 formation, but also 
a lower SFR, resulting in less efficient H 2 destruction). 
Recipes that do not account for the effect of a vary¬ 
ing UV radiation field (GKFUV and KMT) predict that 
H 2 formation becomes so inefficient in low mass galax¬ 
ies that the stellar mass function actually turns over at 
log(m*/M Q ) ~ 8. Similarly, / s tar(Afh) declines sharply 
at \og(Mh/M.Q) ~ 10. Note that although in the models 
shown, the abundance of very low-mass galaxies is actu¬ 
ally probably too small compared with observations, we 
could probably adjust our recipes for stellar feedback or 
photo-ionization squelching to fix this. However, it does 
appear that even in these models, the excess of low-mass 
galaxies (m s tar ~ 10 9-10 Mg) at intermediate redshift 
(0.5 S z S 4) persists. This indicates that the halo mass 
scale where star formation can become inefficient enough 
to break the self-regulation equilibrium is smaller than 
the one where the discrepancy with current observations 
appears. 

The second point is that the differences between 
models seen in more massive halos are almost solely due 
to the assumed scaling of the star formation relation at 
large gas surface densities. The model with the steep¬ 
est dependence of Esfr on Eh 2 (Big2, with Nsf —> 2 
in dense gas) has the highest number density of massive 
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Figure 22. Stellar mass function evolution with redshift. Symbols show observational estimates as detailed in Fig. 1111 The purple 
solid line shows the results of the fiducial GK model, the cyan dotted line shows the GKFUV model, the dashed lavender line 
shows the GK Bigl model and the long-dashed blue line shows the KMT model (see text and Table El. 


galaxies and the highest values of / s t ar in massive halos. 
The GKFUV model is almost identical to the fiducial GK 
model at high masses. The KMT model is the next high¬ 
est (./Vsf —>• 1.4 in dense gas), then the GK Bigl model 
(TVsf = 1). This is because regardless of the gas partion- 
ing recipe, gas in these galaxies is dense enough that it 
is nearly all molecular. Stellar driven winds cannot effi¬ 
ciently escape these deep potential wells. Therefore there 
is a strong dependence on the gas depletion time (star 
formation efficiency). 

ft is also clear from Fig. [23] that modeling of photo¬ 
ionization squelching will have an extremely degenerate 
effect with that of gas partioning and stellar feedback on 
stellar properties. However, observations of gas content 
should help break these degeneracies. 
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Figure 23. The stellar mass of central galaxies divided by the total mass of their dark mat ter halo, in redshift bi ns from z — 0-6. 
Dark gray solid lines and shaded areas show constraints from abundance matching from iBehroozi et alj d2013h . Models shown 
are as in Fig. I22l In addition, the dot-dashed purple line (shown in the z = 0 panel only) shows the fiducial GK model with 
photo-ionization squelching switched off. 
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